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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_interp.F90 in branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 7971

Last change on this file since 7971 was 7971, checked in by jchanut, 7 years ago

Add zstar coordinate with AGRIF

  • Property svn:keywords set to Id
File size: 47.9 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_a(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)*hur_a(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_a(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)*hur_a(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_a(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)*hur_a(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_a(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)*hur_a(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_a(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)*hvr_a(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_a(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)*hvr_a(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_a(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)*hvr_a(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_a(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)*hvr_a(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      ! Interpolate barotropic fluxes
507      Agrif_SpecialValue=0.
508      Agrif_UseSpecialValue = ln_spc_dyn
509
510      IF (ll_int_cons) THEN ! Conservative interpolation
511         ! orders matters here !!!!!!
512         CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1._wp, procname=interpub2b) ! Time integrated
513         CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1._wp, procname=interpvb2b)
514         bdy_tinterp = 1
515         CALL Agrif_Bc_variable(unb_id ,calledweight=1._wp, procname=interpunb) ! After
516         CALL Agrif_Bc_variable(vnb_id ,calledweight=1._wp, procname=interpvnb)
517         bdy_tinterp = 2
518         CALL Agrif_Bc_variable(unb_id ,calledweight=0._wp, procname=interpunb) ! Before
519         CALL Agrif_Bc_variable(vnb_id ,calledweight=0._wp, procname=interpvnb)         
520      ELSE ! Linear interpolation
521         bdy_tinterp = 0
522         ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 
523         ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 
524         ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 
525         ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 
526         CALL Agrif_Bc_variable(unb_id,calledweight=zt, procname=interpunb)
527         CALL Agrif_Bc_variable(vnb_id,calledweight=zt, procname=interpvnb)
528      ENDIF
529      Agrif_UseSpecialValue = .FALSE.
530      !
531   END SUBROUTINE Agrif_dta_ts
532
533   SUBROUTINE Agrif_ssh( kt )
534      !!----------------------------------------------------------------------
535      !!                  ***  ROUTINE Agrif_ssh  ***
536      !!---------------------------------------------------------------------- 
537      INTEGER, INTENT(in) ::   kt
538      !!
539      INTEGER :: ji, jj
540      !!---------------------------------------------------------------------- 
541
542      IF( Agrif_Root() )   RETURN
543
544      ! Linear interpolation of sea level
545      Agrif_SpecialValue    = 0.e0
546      Agrif_UseSpecialValue = .TRUE.
547      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
548      Agrif_UseSpecialValue = .FALSE.
549
550      IF((nbondi == -1).OR.(nbondi == 2)) THEN
551         DO jj=1,jpj
552            ssha(2,jj) = hbdy_w(jj)
553         END DO
554      ENDIF
555
556      IF((nbondi == 1).OR.(nbondi == 2)) THEN
557         DO jj=1,jpj
558            ssha(nlci-1,jj) = hbdy_e(jj)
559         END DO
560      ENDIF
561
562      IF((nbondj == -1).OR.(nbondj == 2)) THEN
563         DO ji=1,jpi
564            ssha(ji,2) = hbdy_s(ji)
565         END DO
566      ENDIF
567
568      IF((nbondj == 1).OR.(nbondj == 2)) THEN
569         DO ji=1,jpi
570            ssha(ji,nlcj-1) = hbdy_n(ji)
571         END DO
572      ENDIF
573
574   END SUBROUTINE Agrif_ssh
575
576   SUBROUTINE Agrif_ssh_ts( jn )
577      !!----------------------------------------------------------------------
578      !!                  ***  ROUTINE Agrif_ssh_ts  ***
579      !!---------------------------------------------------------------------- 
580      INTEGER, INTENT(in) ::   jn
581      !!
582      INTEGER :: ji,jj
583      !!---------------------------------------------------------------------- 
584
585      IF((nbondi == -1).OR.(nbondi == 2)) THEN
586         DO jj=1,jpj
587            ssha_e(2,jj) = hbdy_w(jj)
588         END DO
589      ENDIF
590
591      IF((nbondi == 1).OR.(nbondi == 2)) THEN
592         DO jj=1,jpj
593            ssha_e(nlci-1,jj) = hbdy_e(jj)
594         END DO
595      ENDIF
596
597      IF((nbondj == -1).OR.(nbondj == 2)) THEN
598         DO ji=1,jpi
599            ssha_e(ji,2) = hbdy_s(ji)
600         END DO
601      ENDIF
602
603      IF((nbondj == 1).OR.(nbondj == 2)) THEN
604         DO ji=1,jpi
605            ssha_e(ji,nlcj-1) = hbdy_n(ji)
606         END DO
607      ENDIF
608
609   END SUBROUTINE Agrif_ssh_ts
610
611# if defined key_zdftke
612   SUBROUTINE Agrif_tke
613      !!----------------------------------------------------------------------
614      !!                  ***  ROUTINE Agrif_tke  ***
615      !!---------------------------------------------------------------------- 
616      REAL(wp) ::   zalpha
617      !
618      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
619      IF( zalpha > 1. )   zalpha = 1.
620     
621      Agrif_SpecialValue    = 0.e0
622      Agrif_UseSpecialValue = .TRUE.
623     
624      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
625             
626      Agrif_UseSpecialValue = .FALSE.
627      !
628   END SUBROUTINE Agrif_tke
629# endif
630
631   SUBROUTINE interptsn(ptab,i1,i2,j1,j2,k1,k2,n1,n2,before,nb,ndir)
632      !!---------------------------------------------
633      !!   *** ROUTINE interptsn ***
634      !!---------------------------------------------
635      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
636      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
637      LOGICAL, INTENT(in) :: before
638      INTEGER, INTENT(in) :: nb , ndir
639      !
640      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
641      INTEGER :: imin, imax, jmin, jmax
642      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
643      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
644      LOGICAL :: western_side, eastern_side,northern_side,southern_side
645
646      IF (before) THEN         
647         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
648      ELSE
649         !
650         western_side  = (nb == 1).AND.(ndir == 1)
651         eastern_side  = (nb == 1).AND.(ndir == 2)
652         southern_side = (nb == 2).AND.(ndir == 1)
653         northern_side = (nb == 2).AND.(ndir == 2)
654         !
655         zrhox = Agrif_Rhox()
656         !
657         zalpha1 = ( zrhox - 1. ) * 0.5
658         zalpha2 = 1. - zalpha1
659         !
660         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
661         zalpha4 = 1. - zalpha3
662         !
663         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
664         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
665         zalpha5 = 1. - zalpha6 - zalpha7
666         !
667         imin = i1
668         imax = i2
669         jmin = j1
670         jmax = j2
671         !
672         ! Remove CORNERS
673         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
674         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
675         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
676         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
677         !
678         IF( eastern_side) THEN
679            DO jn = 1, jpts
680               tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
681               DO jk = 1, jpkm1
682                  DO jj = jmin,jmax
683                     IF( umask(nlci-2,jj,jk) == 0.e0 ) THEN
684                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
685                     ELSE
686                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
687                        IF( un(nlci-2,jj,jk) > 0.e0 ) THEN
688                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
689                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
690                        ENDIF
691                     ENDIF
692                  END DO
693               END DO
694            ENDDO
695         ENDIF
696         !
697         IF( northern_side ) THEN           
698            DO jn = 1, jpts
699               tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
700               DO jk = 1, jpkm1
701                  DO ji = imin,imax
702                     IF( vmask(ji,nlcj-2,jk) == 0.e0 ) THEN
703                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
704                     ELSE
705                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
706                        IF (vn(ji,nlcj-2,jk) > 0.e0 ) THEN
707                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
708                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
709                        ENDIF
710                     ENDIF
711                  END DO
712               END DO
713            ENDDO
714         ENDIF
715         !
716         IF( western_side) THEN           
717            DO jn = 1, jpts
718               tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
719               DO jk = 1, jpkm1
720                  DO jj = jmin,jmax
721                     IF( umask(2,jj,jk) == 0.e0 ) THEN
722                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
723                     ELSE
724                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
725                        IF( un(2,jj,jk) < 0.e0 ) THEN
726                           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)
727                        ENDIF
728                     ENDIF
729                  END DO
730               END DO
731            END DO
732         ENDIF
733         !
734         IF( southern_side ) THEN           
735            DO jn = 1, jpts
736               tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
737               DO jk=1,jpk     
738                  DO ji=imin,imax
739                     IF( vmask(ji,2,jk) == 0.e0 ) THEN
740                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
741                     ELSE
742                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
743                        IF( vn(ji,2,jk) < 0.e0 ) THEN
744                           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)
745                        ENDIF
746                     ENDIF
747                  END DO
748               END DO
749            ENDDO
750         ENDIF
751         !
752         ! Treatment of corners
753         !
754         ! East south
755         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
756            tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
757         ENDIF
758         ! East north
759         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
760            tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
761         ENDIF
762         ! West south
763         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
764            tsa(2,2,:,:) = ptab(2,2,:,:)
765         ENDIF
766         ! West north
767         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
768            tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
769         ENDIF
770         !
771      ENDIF
772      !
773   END SUBROUTINE interptsn
774
775   SUBROUTINE interpsshn(ptab,i1,i2,j1,j2,before,nb,ndir)
776      !!----------------------------------------------------------------------
777      !!                  ***  ROUTINE interpsshn  ***
778      !!---------------------------------------------------------------------- 
779      INTEGER, INTENT(in) :: i1,i2,j1,j2
780      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
781      LOGICAL, INTENT(in) :: before
782      INTEGER, INTENT(in) :: nb , ndir
783      LOGICAL :: western_side, eastern_side,northern_side,southern_side
784      !!---------------------------------------------------------------------- 
785      !
786      IF( before) THEN
787         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
788      ELSE
789         western_side  = (nb == 1).AND.(ndir == 1)
790         eastern_side  = (nb == 1).AND.(ndir == 2)
791         southern_side = (nb == 2).AND.(ndir == 1)
792         northern_side = (nb == 2).AND.(ndir == 2)
793         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
794         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
795         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
796         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
797      ENDIF
798      !
799   END SUBROUTINE interpsshn
800
801   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2, before)
802      !!---------------------------------------------
803      !!   *** ROUTINE interpun ***
804      !!---------------------------------------------   
805      !!
806      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
807      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
808      LOGICAL, INTENT(in) :: before
809      !!
810      INTEGER :: ji,jj,jk
811      REAL(wp) :: zrhoy 
812      !!---------------------------------------------   
813      !
814      IF (before) THEN
815         DO jk=1,jpk
816            DO jj=j1,j2
817               DO ji=i1,i2
818                  ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
819                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u_n(ji,jj,jk)
820               END DO
821            END DO
822         END DO
823      ELSE
824         zrhoy = Agrif_Rhoy()
825         DO jk=1,jpkm1
826            DO jj=j1,j2
827               ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj)))
828               ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u_a(i1:i2,jj,jk)
829            END DO
830         END DO
831      ENDIF
832      !
833   END SUBROUTINE interpun
834
835
836   SUBROUTINE interpun2d(ptab,i1,i2,j1,j2,before)
837      !!---------------------------------------------
838      !!   *** ROUTINE interpun ***
839      !!---------------------------------------------   
840      !
841      INTEGER, INTENT(in) :: i1,i2,j1,j2
842      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
843      LOGICAL, INTENT(in) :: before
844      !
845      INTEGER :: ji,jj
846      REAL(wp) :: ztref
847      REAL(wp) :: zrhoy 
848      !!---------------------------------------------   
849      !
850      ztref = 1.
851
852      IF (before) THEN
853         DO jj=j1,j2
854            DO ji=i1,MIN(i2,nlci-1)
855               ptab(ji,jj) = e2u(ji,jj) * ((gcx(ji+1,jj) - gcx(ji,jj))/e1u(ji,jj)) 
856            END DO
857         END DO
858      ELSE
859         zrhoy = Agrif_Rhoy()
860         DO jj=j1,j2
861            laplacu(i1:i2,jj) = ztref * (ptab(i1:i2,jj)/(zrhoy*e2u(i1:i2,jj))) !*umask(i1:i2,jj,1)
862         END DO
863      ENDIF
864      !
865   END SUBROUTINE interpun2d
866
867
868   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before)
869      !!---------------------------------------------
870      !!   *** ROUTINE interpvn ***
871      !!---------------------------------------------   
872      !
873      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
874      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
875      LOGICAL, INTENT(in) :: before
876      !
877      INTEGER :: ji,jj,jk
878      REAL(wp) :: zrhox 
879      !!---------------------------------------------   
880      !     
881      IF (before) THEN         
882         !interpv entre 1 et k2 et interpv2d en jpkp1
883         DO jk=k1,jpk
884            DO jj=j1,j2
885               DO ji=i1,i2
886                  ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
887                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v_n(ji,jj,jk)
888               END DO
889            END DO
890         END DO
891      ELSE         
892         zrhox= Agrif_Rhox()
893         DO jk=1,jpkm1
894            DO jj=j1,j2
895               va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
896               va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v_a(i1:i2,jj,jk)
897            END DO
898         END DO
899      ENDIF
900      !       
901   END SUBROUTINE interpvn
902
903   SUBROUTINE interpvn2d(ptab,i1,i2,j1,j2,before)
904      !!---------------------------------------------
905      !!   *** ROUTINE interpvn ***
906      !!---------------------------------------------   
907      !
908      INTEGER, INTENT(in) :: i1,i2,j1,j2
909      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
910      LOGICAL, INTENT(in) :: before
911      !
912      INTEGER :: ji,jj
913      REAL(wp) :: zrhox 
914      REAL(wp) :: ztref
915      !!---------------------------------------------   
916      !
917      ztref = 1.   
918      IF (before) THEN 
919         !interpv entre 1 et k2 et interpv2d en jpkp1
920         DO jj=j1,MIN(j2,nlcj-1)
921            DO ji=i1,i2
922               ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1)
923            END DO
924         END DO
925      ELSE           
926         zrhox = Agrif_Rhox()
927         DO ji=i1,i2
928            laplacv(ji,j1:j2) = ztref * (ptab(ji,j1:j2)/(zrhox*e1v(ji,j1:j2)))
929         END DO
930      ENDIF
931      !     
932   END SUBROUTINE interpvn2d
933
934   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir)
935      !!----------------------------------------------------------------------
936      !!                  ***  ROUTINE interpunb  ***
937      !!---------------------------------------------------------------------- 
938      INTEGER, INTENT(in) :: i1,i2,j1,j2
939      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
940      LOGICAL, INTENT(in) :: before
941      INTEGER, INTENT(in) :: nb , ndir
942      !!
943      INTEGER :: ji,jj
944      REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
945      LOGICAL :: western_side, eastern_side,northern_side,southern_side
946      !!---------------------------------------------------------------------- 
947      !
948      IF (before) THEN
949         DO jj=j1,j2
950            DO ji=i1,i2
951               ptab(ji,jj) = un_b(ji,jj) * e2u(ji,jj) * hu(ji,jj) 
952            END DO
953         END DO
954      ELSE
955         western_side  = (nb == 1).AND.(ndir == 1)
956         eastern_side  = (nb == 1).AND.(ndir == 2)
957         southern_side = (nb == 2).AND.(ndir == 1)
958         northern_side = (nb == 2).AND.(ndir == 2)
959         zrhoy = Agrif_Rhoy()
960         zrhot = Agrif_rhot()
961         ! Time indexes bounds for integration
962         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
963         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
964         ! Polynomial interpolation coefficients:
965         IF( bdy_tinterp == 1 ) THEN
966            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
967                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
968         ELSEIF( bdy_tinterp == 2 ) THEN
969            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
970                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
971
972         ELSE
973            ztcoeff = 1
974         ENDIF
975         !   
976         IF(western_side) THEN
977            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
978         ENDIF
979         IF(eastern_side) THEN
980            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
981         ENDIF
982         IF(southern_side) THEN
983            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
984         ENDIF
985         IF(northern_side) THEN
986            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
987         ENDIF
988         !           
989         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
990            IF(western_side) THEN
991               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
992                     &                                  * umask(i1,j1:j2,1)
993            ENDIF
994            IF(eastern_side) THEN
995               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
996                     &                                  * umask(i1,j1:j2,1)
997            ENDIF
998            IF(southern_side) THEN
999               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
1000                     &                                  * umask(i1:i2,j1,1)
1001            ENDIF
1002            IF(northern_side) THEN
1003               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
1004                     &                                  * umask(i1:i2,j1,1)
1005            ENDIF
1006         ENDIF
1007      ENDIF
1008      !
1009   END SUBROUTINE interpunb
1010
1011   SUBROUTINE interpvnb(ptab,i1,i2,j1,j2,before,nb,ndir)
1012      !!----------------------------------------------------------------------
1013      !!                  ***  ROUTINE interpvnb  ***
1014      !!---------------------------------------------------------------------- 
1015      INTEGER, INTENT(in) :: i1,i2,j1,j2
1016      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1017      LOGICAL, INTENT(in) :: before
1018      INTEGER, INTENT(in) :: nb , ndir
1019      !!
1020      INTEGER :: ji,jj
1021      REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff   
1022      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1023      !!---------------------------------------------------------------------- 
1024      !
1025      IF (before) THEN
1026         DO jj=j1,j2
1027            DO ji=i1,i2
1028               ptab(ji,jj) = vn_b(ji,jj) * e1v(ji,jj) * hv(ji,jj) 
1029            END DO
1030         END DO
1031      ELSE
1032         western_side  = (nb == 1).AND.(ndir == 1)
1033         eastern_side  = (nb == 1).AND.(ndir == 2)
1034         southern_side = (nb == 2).AND.(ndir == 1)
1035         northern_side = (nb == 2).AND.(ndir == 2)
1036         zrhox = Agrif_Rhox()
1037         zrhot = Agrif_rhot()
1038         ! Time indexes bounds for integration
1039         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1040         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1041         IF( bdy_tinterp == 1 ) THEN
1042            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1043                  &      - zt0**2._wp * (       zt0 - 1._wp)        )
1044         ELSEIF( bdy_tinterp == 2 ) THEN
1045            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1046                  &      - zt0        * (       zt0 - 1._wp)**2._wp ) 
1047
1048         ELSE
1049            ztcoeff = 1
1050         ENDIF
1051         !
1052         IF(western_side) THEN
1053            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1054         ENDIF
1055         IF(eastern_side) THEN
1056            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1057         ENDIF
1058         IF(southern_side) THEN
1059            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1060         ENDIF
1061         IF(northern_side) THEN
1062            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1063         ENDIF
1064         !           
1065         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1066            IF(western_side) THEN
1067               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1068                     &                                  * vmask(i1,j1:j2,1)
1069            ENDIF
1070            IF(eastern_side) THEN
1071               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1072                     &                                  * vmask(i1,j1:j2,1)
1073            ENDIF
1074            IF(southern_side) THEN
1075               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1076                     &                                  * vmask(i1:i2,j1,1)
1077            ENDIF
1078            IF(northern_side) THEN
1079               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1080                     &                                  * vmask(i1:i2,j1,1)
1081            ENDIF
1082         ENDIF
1083      ENDIF
1084      !
1085   END SUBROUTINE interpvnb
1086
1087   SUBROUTINE interpub2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1088      !!----------------------------------------------------------------------
1089      !!                  ***  ROUTINE interpub2b  ***
1090      !!---------------------------------------------------------------------- 
1091      INTEGER, INTENT(in) :: i1,i2,j1,j2
1092      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1093      LOGICAL, INTENT(in) :: before
1094      INTEGER, INTENT(in) :: nb , ndir
1095      !!
1096      INTEGER :: ji,jj
1097      REAL(wp) :: zrhot, zt0, zt1,zat
1098      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1099      !!---------------------------------------------------------------------- 
1100      IF( before ) THEN
1101         DO jj=j1,j2
1102            DO ji=i1,i2
1103               ptab(ji,jj) = ub2_b(ji,jj) * e2u(ji,jj)
1104            END DO
1105         END DO
1106      ELSE
1107         western_side  = (nb == 1).AND.(ndir == 1)
1108         eastern_side  = (nb == 1).AND.(ndir == 2)
1109         southern_side = (nb == 2).AND.(ndir == 1)
1110         northern_side = (nb == 2).AND.(ndir == 2)
1111         zrhot = Agrif_rhot()
1112         ! Time indexes bounds for integration
1113         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1114         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1115         ! Polynomial interpolation coefficients:
1116         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1117                 &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1118         !
1119         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1120         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1121         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1122         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1123      ENDIF
1124      !
1125   END SUBROUTINE interpub2b
1126
1127   SUBROUTINE interpvb2b(ptab,i1,i2,j1,j2,before,nb,ndir)
1128      !!----------------------------------------------------------------------
1129      !!                  ***  ROUTINE interpvb2b  ***
1130      !!---------------------------------------------------------------------- 
1131      INTEGER, INTENT(in) :: i1,i2,j1,j2
1132      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1133      LOGICAL, INTENT(in) :: before
1134      INTEGER, INTENT(in) :: nb , ndir
1135      !!
1136      INTEGER :: ji,jj
1137      REAL(wp) :: zrhot, zt0, zt1,zat
1138      LOGICAL :: western_side, eastern_side,northern_side,southern_side
1139      !!---------------------------------------------------------------------- 
1140      !
1141      IF( before ) THEN
1142         DO jj=j1,j2
1143            DO ji=i1,i2
1144               ptab(ji,jj) = vb2_b(ji,jj) * e1v(ji,jj)
1145            END DO
1146         END DO
1147      ELSE     
1148         western_side  = (nb == 1).AND.(ndir == 1)
1149         eastern_side  = (nb == 1).AND.(ndir == 2)
1150         southern_side = (nb == 2).AND.(ndir == 1)
1151         northern_side = (nb == 2).AND.(ndir == 2)
1152         zrhot = Agrif_rhot()
1153         ! Time indexes bounds for integration
1154         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1155         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1156         ! Polynomial interpolation coefficients:
1157         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
1158                 &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
1159         !
1160         IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1161         IF(eastern_side ) vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1162         IF(southern_side) vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1163         IF(northern_side) vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1164      ENDIF
1165      !     
1166   END SUBROUTINE interpvb2b
1167
1168   SUBROUTINE interpe3t(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1169      !!----------------------------------------------------------------------
1170      !!                  ***  ROUTINE interpe3t  ***
1171      !!---------------------------------------------------------------------- 
1172      !
1173      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1174      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1175      LOGICAL :: before
1176      INTEGER, INTENT(in) :: nb , ndir
1177      !
1178      INTEGER :: ji, jj, jk
1179      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1180      REAL(wp) :: ztmpmsk     
1181      !!---------------------------------------------------------------------- 
1182      !   
1183      IF (before) THEN
1184         DO jk=k1,k2
1185            DO jj=j1,j2
1186               DO ji=i1,i2
1187                  ptab(ji,jj,jk) = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
1188               END DO
1189            END DO
1190         END DO
1191      ELSE
1192         western_side  = (nb == 1).AND.(ndir == 1)
1193         eastern_side  = (nb == 1).AND.(ndir == 2)
1194         southern_side = (nb == 2).AND.(ndir == 1)
1195         northern_side = (nb == 2).AND.(ndir == 2)
1196
1197         DO jk=k1,k2
1198            DO jj=j1,j2
1199               DO ji=i1,i2
1200                  ! Get velocity mask at boundary edge points:
1201                  IF (western_side)  ztmpmsk = umask(ji    ,jj    ,1)
1202                  IF (eastern_side)  ztmpmsk = umask(nlci-2,jj    ,1)
1203                  IF (northern_side) ztmpmsk = vmask(ji    ,nlcj-2,1)
1204                  IF (southern_side) ztmpmsk = vmask(ji    ,2     ,1)
1205
1206                  IF (ABS(ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk))*ztmpmsk > 1.D-2) THEN
1207                     IF (western_side) THEN
1208                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1209                     ELSEIF (eastern_side) THEN
1210                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1211                     ELSEIF (southern_side) THEN
1212                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1213                     ELSEIF (northern_side) THEN
1214                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1215                     ENDIF
1216                     WRITE(numout,*) '      ptab(ji,jj,jk), fse3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1217                     kindic_agr = kindic_agr + 1
1218                  ENDIF
1219               END DO
1220            END DO
1221         END DO
1222
1223      ENDIF
1224      !
1225   END SUBROUTINE interpe3t
1226
1227   SUBROUTINE interpumsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1228      !!----------------------------------------------------------------------
1229      !!                  ***  ROUTINE interpumsk  ***
1230      !!---------------------------------------------------------------------- 
1231      !
1232      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1233      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1234      LOGICAL :: before
1235      INTEGER, INTENT(in) :: nb , ndir
1236      !
1237      INTEGER :: ji, jj, jk
1238      LOGICAL :: western_side, eastern_side   
1239      !!---------------------------------------------------------------------- 
1240      !   
1241      IF (before) THEN
1242         DO jk=k1,k2
1243            DO jj=j1,j2
1244               DO ji=i1,i2
1245                  ptab(ji,jj,jk) = umask(ji,jj,jk)
1246               END DO
1247            END DO
1248         END DO
1249      ELSE
1250
1251         western_side  = (nb == 1).AND.(ndir == 1)
1252         eastern_side  = (nb == 1).AND.(ndir == 2)
1253         DO jk=k1,k2
1254            DO jj=j1,j2
1255               DO ji=i1,i2
1256                   ! Velocity mask at boundary edge points:
1257                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1258                     IF (western_side) THEN
1259                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1260                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1261                        kindic_agr = kindic_agr + 1
1262                     ELSEIF (eastern_side) THEN
1263                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1264                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1265                        kindic_agr = kindic_agr + 1
1266                     ENDIF
1267                  ENDIF
1268               END DO
1269            END DO
1270         END DO
1271
1272      ENDIF
1273      !
1274   END SUBROUTINE interpumsk
1275
1276   SUBROUTINE interpvmsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
1277      !!----------------------------------------------------------------------
1278      !!                  ***  ROUTINE interpvmsk  ***
1279      !!---------------------------------------------------------------------- 
1280      !
1281      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1282      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1283      LOGICAL :: before
1284      INTEGER, INTENT(in) :: nb , ndir
1285      !
1286      INTEGER :: ji, jj, jk
1287      LOGICAL :: northern_side, southern_side     
1288      !!---------------------------------------------------------------------- 
1289      !   
1290      IF (before) THEN
1291         DO jk=k1,k2
1292            DO jj=j1,j2
1293               DO ji=i1,i2
1294                  ptab(ji,jj,jk) = vmask(ji,jj,jk)
1295               END DO
1296            END DO
1297         END DO
1298      ELSE
1299
1300         southern_side = (nb == 2).AND.(ndir == 1)
1301         northern_side = (nb == 2).AND.(ndir == 2)
1302         DO jk=k1,k2
1303            DO jj=j1,j2
1304               DO ji=i1,i2
1305                   ! Velocity mask at boundary edge points:
1306                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1307                     IF (southern_side) THEN
1308                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1309                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1310                        kindic_agr = kindic_agr + 1
1311                     ELSEIF (northern_side) THEN
1312                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1313                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1314                        kindic_agr = kindic_agr + 1
1315                     ENDIF
1316                  ENDIF
1317               END DO
1318            END DO
1319         END DO
1320
1321      ENDIF
1322      !
1323   END SUBROUTINE interpvmsk
1324
1325# if defined key_zdftke
1326
1327   SUBROUTINE interpavm(ptab,i1,i2,j1,j2,k1,k2,before)
1328      !!----------------------------------------------------------------------
1329      !!                  ***  ROUTINE interavm  ***
1330      !!---------------------------------------------------------------------- 
1331      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
1332      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1333      LOGICAL, INTENT(in) :: before
1334      !!---------------------------------------------------------------------- 
1335      !     
1336      IF( before) THEN
1337         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1338      ELSE
1339         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1340      ENDIF
1341      !
1342   END SUBROUTINE interpavm
1343
1344# endif /* key_zdftke */
1345
1346#else
1347   !!----------------------------------------------------------------------
1348   !!   Empty module                                          no AGRIF zoom
1349   !!----------------------------------------------------------------------
1350CONTAINS
1351   SUBROUTINE Agrif_OPA_Interp_empty
1352      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1353   END SUBROUTINE Agrif_OPA_Interp_empty
1354#endif
1355
1356   !!======================================================================
1357END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.