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_oce_interp.F90 in NEMO/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_oce_interp.F90 @ 10165

Last change on this file since 10165 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

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