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

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 6 years ago

Remove svn keywords

File size: 47.8 KB
RevLine 
[636]1MODULE agrif_opa_interp
[1605]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)
[6204]9   !!            3.6  !  2014-09  (R. Benshila)
[1605]10   !!----------------------------------------------------------------------
[2528]11#if defined key_agrif && ! defined key_offline
[1605]12   !!----------------------------------------------------------------------
13   !!   'key_agrif'                                              AGRIF zoom
[2528]14   !!   NOT 'key_offline'                               NO off-line tracers
[1605]15   !!----------------------------------------------------------------------
16   !!   Agrif_tra     :
17   !!   Agrif_dyn     :
18   !!   interpu       :
19   !!   interpv       :
20   !!----------------------------------------------------------------------
[636]21   USE par_oce
22   USE oce
23   USE dom_oce     
24   USE sol_oce
[782]25   USE agrif_oce
[1605]26   USE phycst
27   USE in_out_manager
[2715]28   USE agrif_opa_sponge
29   USE lib_mpp
[4292]30   USE wrk_nemo
[4486]31   USE dynspg_oce
[6204]32   USE zdf_oce
33 
[636]34   IMPLICIT NONE
35   PRIVATE
[4292]36
[6204]37   INTEGER :: bdy_tinterp = 0
38
[4486]39   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
[6204]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
[390]47
[1605]48#  include "domzgr_substitute.h90" 
49#  include "vectopt_loop_substitute.h90"
[1156]50   !!----------------------------------------------------------------------
[6204]51   !! NEMO/NST 3.6 , NEMO Consortium (2010)
[1156]52   !! $Id$
[2528]53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1156]54   !!----------------------------------------------------------------------
55
[6204]56CONTAINS
57
[782]58   SUBROUTINE Agrif_tra
[1605]59      !!----------------------------------------------------------------------
[6204]60      !!                  ***  ROUTINE Agrif_tra  ***
[1605]61      !!----------------------------------------------------------------------
[636]62      !
[1605]63      IF( Agrif_Root() )   RETURN
[390]64
[1605]65      Agrif_SpecialValue    = 0.e0
[636]66      Agrif_UseSpecialValue = .TRUE.
[390]67
[6204]68      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
[636]69      Agrif_UseSpecialValue = .FALSE.
[1605]70      !
[636]71   END SUBROUTINE Agrif_tra
72
[1605]73
[636]74   SUBROUTINE Agrif_dyn( kt )
[1605]75      !!----------------------------------------------------------------------
76      !!                  ***  ROUTINE Agrif_DYN  ***
77      !!---------------------------------------------------------------------- 
[2715]78      !!
[1605]79      INTEGER, INTENT(in) ::   kt
80      !!
[6204]81      INTEGER :: ji,jj,jk, j1,j2, i1,i2
[636]82      REAL(wp) :: timeref
[390]83      REAL(wp) :: z2dt, znugdt
[4292]84      REAL(wp) :: zrhox, zrhoy
[6204]85      REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1
[1605]86      !!---------------------------------------------------------------------- 
[390]87
[1605]88      IF( Agrif_Root() )   RETURN
[390]89
[6204]90      CALL wrk_alloc( jpi, jpj, spgv1, spgu1 )
[2715]91
[6204]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
[636]105      zrhox = Agrif_Rhox()
[4292]106      zrhoy = Agrif_Rhoy()
[390]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
[1605]114      znugdt =  grav * z2dt   
[390]115
[6204]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
[782]125
[390]126
[636]127      IF((nbondi == -1).OR.(nbondi == 2)) THEN
[4292]128#if defined key_dynspg_flt
[636]129         DO jk=1,jpkm1
[6204]130            DO jj=j1,j2
[636]131               ua(2,jj,jk) = (ua(2,jj,jk) - z2dt * znugdt * laplacu(2,jj))*umask(2,jj,jk)
132            END DO
133         END DO
[390]134
[636]135         spgu(2,:)=0.
[390]136
[636]137         DO jk=1,jpkm1
138            DO jj=1,jpj
[6204]139               spgu(2,jj)=spgu(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
[636]140            END DO
141         END DO
[390]142
[636]143         DO jj=1,jpj
144            IF (umask(2,jj,1).NE.0.) THEN
[6204]145               spgu(2,jj)=spgu(2,jj)/hu(2,jj)
[636]146            ENDIF
147         END DO
[4292]148#else
149         spgu(2,:) = ua_b(2,:)
150#endif
[390]151
[636]152         DO jk=1,jpkm1
[6204]153            DO jj=j1,j2
[636]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
[390]158
[636]159         spgu1(2,:)=0.
[390]160
[636]161         DO jk=1,jpkm1
162            DO jj=1,jpj
[6204]163               spgu1(2,jj)=spgu1(2,jj)+fse3u(2,jj,jk)*ua(2,jj,jk)
[636]164            END DO
165         END DO
[390]166
[636]167         DO jj=1,jpj
168            IF (umask(2,jj,1).NE.0.) THEN
[6204]169               spgu1(2,jj)=spgu1(2,jj)/hu(2,jj)
[636]170            ENDIF
171         END DO
[390]172
[636]173         DO jk=1,jpkm1
[6204]174            DO jj=j1,j2
[636]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
[390]178
[4486]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
[636]197      ENDIF
[390]198
[636]199      IF((nbondi == 1).OR.(nbondi == 2)) THEN
[4292]200#if defined key_dynspg_flt
[636]201         DO jk=1,jpkm1
[6204]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)
[636]204            END DO
205         END DO
[6204]206         spgu(nlci-2,:)=0.
[636]207         DO jk=1,jpkm1
208            DO jj=1,jpj
[6204]209               spgu(nlci-2,jj)=spgu(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)
210            ENDDO
211         ENDDO
[636]212         DO jj=1,jpj
213            IF (umask(nlci-2,jj,1).NE.0.) THEN
[6204]214               spgu(nlci-2,jj)=spgu(nlci-2,jj)/hu(nlci-2,jj)
[636]215            ENDIF
216         END DO
[4292]217#else
218         spgu(nlci-2,:) = ua_b(nlci-2,:)
219#endif
[636]220         DO jk=1,jpkm1
[6204]221            DO jj=j1,j2
[636]222               ua(nlci-2,jj,jk) = 0.25*(ua(nlci-3,jj,jk)+2.*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk))
[390]223
[636]224               ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk)
[390]225
[636]226            END DO
227         END DO
228         spgu1(nlci-2,:)=0.
229         DO jk=1,jpkm1
230            DO jj=1,jpj
[6204]231               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)+fse3u(nlci-2,jj,jk)*ua(nlci-2,jj,jk)*umask(nlci-2,jj,jk)
[636]232            END DO
233         END DO
234         DO jj=1,jpj
235            IF (umask(nlci-2,jj,1).NE.0.) THEN
[6204]236               spgu1(nlci-2,jj)=spgu1(nlci-2,jj)/hu(nlci-2,jj)
[636]237            ENDIF
238         END DO
239         DO jk=1,jpkm1
[6204]240            DO jj=j1,j2
[636]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
[390]244
[4486]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
[636]265      ENDIF
[390]266
[636]267      IF((nbondj == -1).OR.(nbondj == 2)) THEN
[390]268
[4292]269#if defined key_dynspg_flt
[636]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
[390]275
[636]276         spgv(:,2)=0.
[390]277
[636]278         DO jk=1,jpkm1
279            DO ji=1,jpi
[6204]280               spgv(ji,2)=spgv(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)
[636]281            END DO
282         END DO
[390]283
[636]284         DO ji=1,jpi
285            IF (vmask(ji,2,1).NE.0.) THEN
[6204]286               spgv(ji,2)=spgv(ji,2)/hv(ji,2)
[636]287            ENDIF
288         END DO
[4292]289#else
290         spgv(:,2)=va_b(:,2)
291#endif
[390]292
[636]293         DO jk=1,jpkm1
[6204]294            DO ji=i1,i2
[636]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
[390]299
[636]300         spgv1(:,2)=0.
[390]301
[636]302         DO jk=1,jpkm1
303            DO ji=1,jpi
[6204]304               spgv1(ji,2)=spgv1(ji,2)+fse3v(ji,2,jk)*va(ji,2,jk)*vmask(ji,2,jk)
[636]305            END DO
306         END DO
[390]307
[636]308         DO ji=1,jpi
309            IF (vmask(ji,2,1).NE.0.) THEN
[6204]310               spgv1(ji,2)=spgv1(ji,2)/hv(ji,2)
[636]311            ENDIF
312         END DO
[390]313
[636]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
[390]319
[4486]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
[636]339      ENDIF
[390]340
[636]341      IF((nbondj == 1).OR.(nbondj == 2)) THEN
[390]342
[4292]343#if defined key_dynspg_flt
[636]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
[390]349
[6204]350
[636]351         spgv(:,nlcj-2)=0.
[390]352
[636]353         DO jk=1,jpkm1
354            DO ji=1,jpi
[6204]355               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
[636]356            END DO
357         END DO
[390]358
[636]359         DO ji=1,jpi
360            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
[6204]361               spgv(ji,nlcj-2)=spgv(ji,nlcj-2)/hv(ji,nlcj-2)
[636]362            ENDIF
363         END DO
[6204]364
[4292]365#else
366         spgv(:,nlcj-2)=va_b(:,nlcj-2)
367#endif
[390]368
[636]369         DO jk=1,jpkm1
[6204]370            DO ji=i1,i2
[636]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
[390]375
[636]376         spgv1(:,nlcj-2)=0.
[390]377
[636]378         DO jk=1,jpkm1
379            DO ji=1,jpi
[6204]380               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)+fse3v(ji,nlcj-2,jk)*va(ji,nlcj-2,jk)
[636]381            END DO
382         END DO
[390]383
[636]384         DO ji=1,jpi
385            IF (vmask(ji,nlcj-2,1).NE.0.) THEN
[6204]386               spgv1(ji,nlcj-2)=spgv1(ji,nlcj-2)/hv(ji,nlcj-2)
[636]387            ENDIF
388         END DO
[390]389
[636]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
[390]395
[4486]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
[636]416      ENDIF
[2715]417      !
[6204]418      CALL wrk_dealloc( jpi, jpj, spgv1, spgu1 )
[2715]419      !
[636]420   END SUBROUTINE Agrif_dyn
[390]421
[4486]422   SUBROUTINE Agrif_dyn_ts( jn )
[4292]423      !!----------------------------------------------------------------------
424      !!                  ***  ROUTINE Agrif_dyn_ts  ***
425      !!---------------------------------------------------------------------- 
426      !!
[4486]427      INTEGER, INTENT(in) ::   jn
[4292]428      !!
429      INTEGER :: ji, jj
[4486]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)
[6204]437            ! Specified fluxes:
[4486]438            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
[6204]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)) )
[4486]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)
[6204]448            ! Specified fluxes:
[4486]449            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
[6204]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)) )
[4486]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)
[6204]459            ! Specified fluxes:
[4486]460            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
[6204]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)) )
[4486]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)
[6204]470            ! Specified fluxes:
[4486]471            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
[6204]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)) )
[4486]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
[6204]489      REAL(wp) :: zrhot, zt
[4292]490      !!---------------------------------------------------------------------- 
[1605]491
[4292]492      IF( Agrif_Root() )   RETURN
493
[4486]494      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in
[6204]495      ! the forward case only
[4486]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
[4292]505
[4486]506      ! Linear interpolation of sea level
507      Agrif_SpecialValue    = 0.e0
508      Agrif_UseSpecialValue = .TRUE.
[6204]509      CALL Agrif_Bc_variable(sshn_id,calledweight=zt, procname=interpsshn )
[4486]510      Agrif_UseSpecialValue = .FALSE.
[4292]511
[4486]512      ! Interpolate barotropic fluxes
513      Agrif_SpecialValue=0.
514      Agrif_UseSpecialValue = ln_spc_dyn
[4292]515
[4486]516      IF (ll_int_cons) THEN ! Conservative interpolation
[6204]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)         
[4486]526      ELSE ! Linear interpolation
[6204]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)
[4486]534      ENDIF
535      Agrif_UseSpecialValue = .FALSE.
[6204]536      !
[4486]537   END SUBROUTINE Agrif_dta_ts
538
[2486]539   SUBROUTINE Agrif_ssh( kt )
540      !!----------------------------------------------------------------------
[2528]541      !!                  ***  ROUTINE Agrif_DYN  ***
[2486]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,:)
[6204]556         sshn(nlci-1,:)=sshn(nlci-2,:)
[2486]557      ENDIF
558
559      IF((nbondj == -1).OR.(nbondj == 2)) THEN
[4292]560         ssha(:,2)=ssha(:,3)
561         sshn(:,2)=sshn(:,3)
[2486]562      ENDIF
563
564      IF((nbondj == 1).OR.(nbondj == 2)) THEN
565         ssha(:,nlcj-1)=ssha(:,nlcj-2)
[6204]566         sshn(:,nlcj-1)=sshn(:,nlcj-2)
[2486]567      ENDIF
568
569   END SUBROUTINE Agrif_ssh
570
[4486]571   SUBROUTINE Agrif_ssh_ts( jn )
[4292]572      !!----------------------------------------------------------------------
573      !!                  ***  ROUTINE Agrif_ssh_ts  ***
574      !!---------------------------------------------------------------------- 
[4486]575      INTEGER, INTENT(in) ::   jn
[4292]576      !!
[4486]577      INTEGER :: ji,jj
[4292]578      !!---------------------------------------------------------------------- 
[2486]579
[4292]580      IF((nbondi == -1).OR.(nbondi == 2)) THEN
[4486]581         DO jj=1,jpj
582            ssha_e(2,jj) = hbdy_w(jj)
583         END DO
[4292]584      ENDIF
585
586      IF((nbondi == 1).OR.(nbondi == 2)) THEN
[4486]587         DO jj=1,jpj
588            ssha_e(nlci-1,jj) = hbdy_e(jj)
589         END DO
[4292]590      ENDIF
591
592      IF((nbondj == -1).OR.(nbondj == 2)) THEN
[4486]593         DO ji=1,jpi
594            ssha_e(ji,2) = hbdy_s(ji)
595         END DO
[4292]596      ENDIF
597
598      IF((nbondj == 1).OR.(nbondj == 2)) THEN
[4486]599         DO ji=1,jpi
600            ssha_e(ji,nlcj-1) = hbdy_n(ji)
601         END DO
[4292]602      ENDIF
603
604   END SUBROUTINE Agrif_ssh_ts
605
[6204]606# if defined key_zdftke
607   SUBROUTINE Agrif_tke
[4292]608      !!----------------------------------------------------------------------
[6204]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      !!----------------------------------------------------------------------
[4292]772      !!                  ***  ROUTINE interpsshn  ***
773      !!---------------------------------------------------------------------- 
774      INTEGER, INTENT(in) :: i1,i2,j1,j2
[6204]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      !!---------------------------------------------   
[4292]800      !!
[6204]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      !
[4292]840      INTEGER :: ji,jj
[6204]841      REAL(wp) :: ztref
842      REAL(wp) :: zrhoy 
843      !!---------------------------------------------   
844      !
845      ztref = 1.
[4292]846
[6204]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
[4292]861
862
[6204]863   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before)
864      !!---------------------------------------------
865      !!   *** ROUTINE interpvn ***
866      !!---------------------------------------------   
867      !
[636]868      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
[6204]869      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
870      LOGICAL, INTENT(in) :: before
871      !
[636]872      INTEGER :: ji,jj,jk
[6204]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
[636]897
[6204]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)
[636]916            DO ji=i1,i2
[6204]917               ptab(ji,jj) = e1v(ji,jj) * ((gcx(ji,jj+1) - gcx(ji,jj))/e2v(ji,jj)) * vmask(ji,jj,1)
[636]918            END DO
919         END DO
[6204]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
[390]928
[6204]929   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir)
[1605]930      !!----------------------------------------------------------------------
[6204]931      !!                  ***  ROUTINE interpunb  ***
[1605]932      !!---------------------------------------------------------------------- 
[636]933      INTEGER, INTENT(in) :: i1,i2,j1,j2
[6204]934      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
935      LOGICAL, INTENT(in) :: before
936      INTEGER, INTENT(in) :: nb , ndir
[1605]937      !!
[636]938      INTEGER :: ji,jj
[6204]939      REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
940      LOGICAL :: western_side, eastern_side,northern_side,southern_side
[1605]941      !!---------------------------------------------------------------------- 
[6204]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
[636]948         END DO
[6204]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 ) 
[636]966
[6204]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
[636]1005
[6204]1006   SUBROUTINE interpvnb(ptab,i1,i2,j1,j2,before,nb,ndir)
[1605]1007      !!----------------------------------------------------------------------
[6204]1008      !!                  ***  ROUTINE interpvnb  ***
[1605]1009      !!---------------------------------------------------------------------- 
[6204]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
[1605]1014      !!
[6204]1015      INTEGER :: ji,jj
1016      REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff   
1017      LOGICAL :: western_side, eastern_side,northern_side,southern_side
[1605]1018      !!---------------------------------------------------------------------- 
[6204]1019      !
1020      IF (before) THEN
[636]1021         DO jj=j1,j2
1022            DO ji=i1,i2
[6204]1023               ptab(ji,jj) = vn_b(ji,jj) * e1v(ji,jj) * hv(ji,jj) 
[636]1024            END DO
1025         END DO
[6204]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 ) 
[390]1042
[6204]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
[390]1081
[6204]1082   SUBROUTINE interpub2b(ptab,i1,i2,j1,j2,before,nb,ndir)
[1605]1083      !!----------------------------------------------------------------------
[6204]1084      !!                  ***  ROUTINE interpub2b  ***
[1605]1085      !!---------------------------------------------------------------------- 
[636]1086      INTEGER, INTENT(in) :: i1,i2,j1,j2
[6204]1087      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1088      LOGICAL, INTENT(in) :: before
1089      INTEGER, INTENT(in) :: nb , ndir
[1605]1090      !!
[636]1091      INTEGER :: ji,jj
[6204]1092      REAL(wp) :: zrhot, zt0, zt1,zat
1093      LOGICAL :: western_side, eastern_side,northern_side,southern_side
[1605]1094      !!---------------------------------------------------------------------- 
[6204]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
[636]1100         END DO
[6204]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
[636]1121
[6204]1122   SUBROUTINE interpvb2b(ptab,i1,i2,j1,j2,before,nb,ndir)
[4292]1123      !!----------------------------------------------------------------------
[6204]1124      !!                  ***  ROUTINE interpvb2b  ***
[4292]1125      !!---------------------------------------------------------------------- 
1126      INTEGER, INTENT(in) :: i1,i2,j1,j2
[6204]1127      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1128      LOGICAL, INTENT(in) :: before
1129      INTEGER, INTENT(in) :: nb , ndir
[4292]1130      !!
[4486]1131      INTEGER :: ji,jj
[6204]1132      REAL(wp) :: zrhot, zt0, zt1,zat
1133      LOGICAL :: western_side, eastern_side,northern_side,southern_side
[4292]1134      !!---------------------------------------------------------------------- 
[6204]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
[4292]1162
[6204]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
[4292]1185         END DO
[6204]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)
[4292]1191
[6204]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)
[4292]1200
[6204]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)
[4292]1223      !!----------------------------------------------------------------------
[6204]1224      !!                  ***  ROUTINE interpumsk  ***
[4292]1225      !!---------------------------------------------------------------------- 
[6204]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   
[4292]1234      !!---------------------------------------------------------------------- 
[6204]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
[4292]1245
[6204]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
[4292]1265         END DO
1266
[6204]1267      ENDIF
1268      !
1269   END SUBROUTINE interpumsk
[4292]1270
[6204]1271   SUBROUTINE interpvmsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
[4486]1272      !!----------------------------------------------------------------------
[6204]1273      !!                  ***  ROUTINE interpvmsk  ***
[4486]1274      !!---------------------------------------------------------------------- 
[6204]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     
[4486]1283      !!---------------------------------------------------------------------- 
[6204]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
[4486]1294
[6204]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
[4486]1314         END DO
1315
[6204]1316      ENDIF
1317      !
1318   END SUBROUTINE interpvmsk
[4486]1319
[6204]1320# if defined key_zdftke
1321
1322   SUBROUTINE interpavm(ptab,i1,i2,j1,j2,k1,k2,before)
[4486]1323      !!----------------------------------------------------------------------
[6204]1324      !!                  ***  ROUTINE interavm  ***
[4486]1325      !!---------------------------------------------------------------------- 
[6204]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
[4486]1329      !!---------------------------------------------------------------------- 
[6204]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
[4486]1338
[6204]1339# endif /* key_zdftke */
[4486]1340
[390]1341#else
[1605]1342   !!----------------------------------------------------------------------
1343   !!   Empty module                                          no AGRIF zoom
1344   !!----------------------------------------------------------------------
[636]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
[390]1349#endif
[1605]1350
1351   !!======================================================================
[636]1352END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.