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.
Changeset 12377 for NEMO/trunk/src/NST/agrif_oce_interp.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/NST/agrif_oce_interp.F90

    r10068 r12377  
    3333   USE agrif_oce_sponge 
    3434   USE lib_mpp 
     35   USE vremap 
    3536  
    3637   IMPLICIT NONE 
    3738   PRIVATE 
    3839 
    39    PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts 
     40   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts 
    4041   PUBLIC   Agrif_tra, Agrif_avm 
    4142   PUBLIC   interpun , interpvn 
    4243   PUBLIC   interptsn, interpsshn, interpavm 
    4344   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b 
    44    PUBLIC   interpe3t, interpumsk, interpvmsk 
    45  
     45   PUBLIC   interpe3t 
     46#if defined key_vertical 
     47   PUBLIC   interpht0, interpmbkt 
     48# endif 
    4649   INTEGER ::   bdy_tinterp = 0 
    4750 
    48 #  include "vectopt_loop_substitute.h90" 
    4951   !!---------------------------------------------------------------------- 
    5052   !! NEMO/NST 4.0 , NEMO Consortium (2018) 
    … …  
    7880      ! 
    7981      INTEGER ::   ji, jj, jk       ! dummy loop indices 
    80       INTEGER ::   j1, j2, i1, i2 
    8182      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2 
    8283      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb 
    … …  
    9394      Agrif_UseSpecialValue = .FALSE. 
    9495      ! 
    95       ! prevent smoothing in ghost cells 
    96       i1 =  1   ;   i2 = nlci 
    97       j1 =  1   ;   j2 = nlcj 
    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 
    102  
    10396      ! --- West --- ! 
    104       IF( nbondi == -1 .OR. nbondi == 2 ) THEN 
    105          ibdy1 = 2 
    106          ibdy2 = 1+nbghostcells  
    107          ! 
    108          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    109             ua_b(ibdy1:ibdy2,:) = 0._wp 
     97      ibdy1 = 2 
     98      ibdy2 = 1+nbghostcells  
     99      ! 
     100      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     101         DO ji = mi0(ibdy1), mi1(ibdy2) 
     102            uu_b(ji,:,Krhs_a) = 0._wp 
     103 
    110104            DO jk = 1, jpkm1 
    111105               DO jj = 1, jpj 
    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) 
    114                END DO 
    115             END DO 
     106                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
     107               END DO 
     108            END DO 
     109 
    116110            DO jj = 1, jpj 
    117                ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
    118             END DO 
    119          ENDIF 
    120          ! 
    121          IF( .NOT.lk_agrif_clp ) THEN 
    122             DO jk=1,jpkm1              ! Smooth 
    123                DO jj=j1,j2 
    124                   ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk)) 
    125                END DO 
    126             END DO 
    127          ENDIF 
    128          ! 
    129          zub(ibdy1:ibdy2,:) = 0._wp    ! Correct transport 
     111               uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     112            END DO 
     113         END DO 
     114      ENDIF 
     115      ! 
     116      DO ji = mi0(ibdy1), mi1(ibdy2) 
     117         zub(ji,:) = 0._wp    ! Correct transport 
    130118         DO jk = 1, jpkm1 
    131119            DO jj = 1, jpj 
    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) 
     120               zub(ji,jj) = zub(ji,jj) &  
     121                  & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a)*umask(ji,jj,jk) 
    134122            END DO 
    135123         END DO 
    136124         DO jj=1,jpj 
    137             zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
     125            zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    138126         END DO 
    139127             
    140128         DO jk = 1, jpkm1 
    141129            DO jj = 1, jpj 
    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) 
    144             END DO 
    145          END DO 
     130               uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) + uu_b(ji,jj,Krhs_a)-zub(ji,jj)) * umask(ji,jj,jk) 
     131            END DO 
     132         END DO 
     133      END DO 
    146134             
    147          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    148             zvb(ibdy1:ibdy2,:) = 0._wp 
     135      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     136         DO ji = mi0(ibdy1), mi1(ibdy2) 
     137            zvb(ji,:) = 0._wp 
    149138            DO jk = 1, jpkm1 
    150139               DO jj = 1, jpj 
    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) 
     140                  zvb(ji,jj) = zvb(ji,jj) + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    153141               END DO 
    154142            END DO 
    155143            DO jj = 1, jpj 
    156                zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj) 
     144               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    157145            END DO 
    158146            DO jk = 1, jpkm1 
    159147               DO jj = 1, jpj 
    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) 
    162                END DO 
    163             END DO 
    164          ENDIF 
    165          ! 
    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  
     148                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) + vv_b(ji,jj,Krhs_a)-zvb(ji,jj))*vmask(ji,jj,jk) 
     149               END DO 
     150            END DO 
     151         END DO 
    172152      ENDIF 
    173153 
    174154      ! --- East --- ! 
    175       IF( nbondi ==  1 .OR. nbondi == 2 ) THEN 
    176          ibdy1 = nlci-1-nbghostcells 
    177          ibdy2 = nlci-2  
    178          ! 
    179          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    180             ua_b(ibdy1:ibdy2,:) = 0._wp 
     155      ibdy1 = jpiglo-1-nbghostcells 
     156      ibdy2 = jpiglo-2  
     157      ! 
     158      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     159         DO ji = mi0(ibdy1), mi1(ibdy2) 
     160            uu_b(ji,:,Krhs_a) = 0._wp 
    181161            DO jk = 1, jpkm1 
    182162               DO jj = 1, jpj 
    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) 
     163                  uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) &  
     164                      & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    185165               END DO 
    186166            END DO 
    187167            DO jj = 1, jpj 
    188                ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
    189             END DO 
    190          ENDIF 
    191          ! 
    192          IF( .NOT.lk_agrif_clp ) THEN 
    193             DO jk=1,jpkm1              ! Smooth 
    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)) 
    196                END DO 
    197             END DO 
    198          ENDIF 
    199          ! 
    200          zub(ibdy1:ibdy2,:) = 0._wp    ! Correct transport 
     168               uu_b(ji,jj,Krhs_a) = uu_b(ji,jj,Krhs_a) * r1_hu(ji,jj,Krhs_a) 
     169            END DO 
     170         END DO 
     171      ENDIF 
     172      ! 
     173      DO ji = mi0(ibdy1), mi1(ibdy2) 
     174         zub(ji,:) = 0._wp    ! Correct transport 
    201175         DO jk = 1, jpkm1 
    202176            DO jj = 1, jpj 
    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) 
     177               zub(ji,jj) = zub(ji,jj) &  
     178                  & + e3u(ji,jj,jk,Krhs_a)  * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    205179            END DO 
    206180         END DO 
    207181         DO jj=1,jpj 
    208             zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
     182            zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    209183         END DO 
    210184             
    211185         DO jk = 1, jpkm1 
    212186            DO jj = 1, jpj 
    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) 
    215             END DO 
    216          END DO 
     187               uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) &  
     188                 & + uu_b(ji,jj,Krhs_a)-zub(ji,jj))*umask(ji,jj,jk) 
     189            END DO 
     190         END DO 
     191      END DO 
    217192             
    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 
     193      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     194         ibdy1 = jpiglo-nbghostcells 
     195         ibdy2 = jpiglo-1  
     196         DO ji = mi0(ibdy1), mi1(ibdy2) 
     197            zvb(ji,:) = 0._wp 
    222198            DO jk = 1, jpkm1 
    223199               DO jj = 1, jpj 
    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) 
     200                  zvb(ji,jj) = zvb(ji,jj) & 
     201                     & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    226202               END DO 
    227203            END DO 
    228204            DO jj = 1, jpj 
    229                zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj) 
     205               zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    230206            END DO 
    231207            DO jk = 1, jpkm1 
    232208               DO jj = 1, jpj 
    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) 
    235                END DO 
    236             END DO 
    237          ENDIF 
    238          ! 
    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  
     209                  vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) &  
     210                      & + vv_b(ji,jj,Krhs_a)-zvb(ji,jj)) * vmask(ji,jj,jk) 
     211               END DO 
     212            END DO 
     213         END DO 
    245214      ENDIF 
    246215 
    247216      ! --- South --- ! 
    248       IF( nbondj == -1 .OR. nbondj == 2 ) THEN 
    249          jbdy1 = 2 
    250          jbdy2 = 1+nbghostcells  
    251          ! 
    252          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    253             va_b(:,jbdy1:jbdy2) = 0._wp 
     217      jbdy1 = 2 
     218      jbdy2 = 1+nbghostcells  
     219      ! 
     220      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     221         DO jj = mj0(jbdy1), mj1(jbdy2) 
     222            vv_b(:,jj,Krhs_a) = 0._wp 
    254223            DO jk = 1, jpkm1 
    255224               DO ji = 1, jpi 
    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) 
     225                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) &  
     226                      & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    258227               END DO 
    259228            END DO 
    260229            DO ji=1,jpi 
    261                va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
    262             END DO 
    263          ENDIF 
    264          ! 
    265          IF ( .NOT.lk_agrif_clp ) THEN 
    266             DO jk = 1, jpkm1           ! Smooth 
    267                DO ji = i1, i2 
    268                   va(ji,jbdy2,jk) = 0.25_wp*(va(ji,jbdy2-1,jk)+2._wp*va(ji,jbdy2,jk)+va(ji,jbdy2+1,jk)) 
    269                END DO 
    270             END DO 
    271          ENDIF 
    272          ! 
    273          zvb(:,jbdy1:jbdy2) = 0._wp    ! Correct transport 
     230               vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a)      
     231            END DO 
     232         END DO 
     233      ENDIF 
     234      ! 
     235      DO jj = mj0(jbdy1), mj1(jbdy2) 
     236         zvb(:,jj) = 0._wp    ! Correct transport 
    274237         DO jk=1,jpkm1 
    275238            DO ji=1,jpi 
    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) 
     239               zvb(ji,jj) = zvb(ji,jj) &  
     240                  & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    278241            END DO 
    279242         END DO 
    280243         DO ji = 1, jpi 
    281             zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
     244            zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    282245         END DO 
    283246 
    284247         DO jk = 1, jpkm1 
    285248            DO ji = 1, jpi 
    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) 
    288             END DO 
    289          END DO 
     249               vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) &  
     250                 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 
     251            END DO 
     252         END DO 
     253      END DO 
    290254             
    291          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    292             zub(:,jbdy1:jbdy2) = 0._wp 
     255      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     256         DO jj = mj0(jbdy1), mj1(jbdy2) 
     257            zub(:,jj) = 0._wp 
    293258            DO jk = 1, jpkm1 
    294259               DO ji = 1, jpi 
    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) 
     260                  zub(ji,jj) = zub(ji,jj) &  
     261                     & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    297262               END DO 
    298263            END DO 
    299264            DO ji = 1, jpi 
    300                zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2) 
     265               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    301266            END DO 
    302267                
    303268            DO jk = 1, jpkm1 
    304269               DO ji = 1, jpi 
    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) 
    307                END DO 
    308             END DO 
    309          ENDIF 
    310          ! 
    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  
     270                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) &  
     271                    & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 
     272               END DO 
     273            END DO 
     274         END DO 
    317275      ENDIF 
    318276 
    319277      ! --- North --- ! 
    320       IF( nbondj ==  1 .OR. nbondj == 2 ) THEN 
    321          jbdy1 = nlcj-1-nbghostcells 
    322          jbdy2 = nlcj-2  
    323          ! 
    324          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
    325             va_b(:,jbdy1:jbdy2) = 0._wp 
     278      jbdy1 = jpjglo-1-nbghostcells 
     279      jbdy2 = jpjglo-2  
     280      ! 
     281      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     282         DO jj = mj0(jbdy1), mj1(jbdy2) 
     283            vv_b(:,jj,Krhs_a) = 0._wp 
    326284            DO jk = 1, jpkm1 
    327285               DO ji = 1, jpi 
    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) 
     286                  vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) &  
     287                      & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    330288               END DO 
    331289            END DO 
    332290            DO ji=1,jpi 
    333                va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
    334             END DO 
    335          ENDIF 
    336          ! 
    337          IF ( .NOT.lk_agrif_clp ) THEN 
    338             DO jk = 1, jpkm1           ! Smooth 
    339                DO ji = i1, i2 
    340                   va(ji,jbdy1,jk) = 0.25_wp*(va(ji,jbdy1-1,jk)+2._wp*va(ji,jbdy1,jk)+va(ji,jbdy1+1,jk)) 
    341                END DO 
    342             END DO 
    343          ENDIF 
    344          ! 
    345          zvb(:,jbdy1:jbdy2) = 0._wp    ! Correct transport 
     291               vv_b(ji,jj,Krhs_a) = vv_b(ji,jj,Krhs_a) * r1_hv(ji,jj,Krhs_a) 
     292            END DO 
     293         END DO 
     294      ENDIF 
     295      ! 
     296      DO jj = mj0(jbdy1), mj1(jbdy2) 
     297         zvb(:,jj) = 0._wp    ! Correct transport 
    346298         DO jk=1,jpkm1 
    347299            DO ji=1,jpi 
    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) 
     300               zvb(ji,jj) = zvb(ji,jj) &  
     301                  & + e3v(ji,jj,jk,Krhs_a) * vv(ji,jj,jk,Krhs_a) * vmask(ji,jj,jk) 
    350302            END DO 
    351303         END DO 
    352304         DO ji = 1, jpi 
    353             zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
     305            zvb(ji,jj) = zvb(ji,jj) * r1_hv(ji,jj,Krhs_a) 
    354306         END DO 
    355307 
    356308         DO jk = 1, jpkm1 
    357309            DO ji = 1, jpi 
    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) 
    360             END DO 
    361          END DO 
     310               vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) &  
     311                 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 
     312            END DO 
     313         END DO 
     314      END DO 
    362315             
    363          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    364             jbdy1 = jbdy1 + 1 
    365             jbdy2 = jbdy2 + 1  
    366             zub(:,jbdy1:jbdy2) = 0._wp 
     316      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     317         jbdy1 = jpjglo-nbghostcells 
     318         jbdy2 = jpjglo-1 
     319         DO jj = mj0(jbdy1), mj1(jbdy2) 
     320            zub(:,jj) = 0._wp 
    367321            DO jk = 1, jpkm1 
    368322               DO ji = 1, jpi 
    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) 
     323                  zub(ji,jj) = zub(ji,jj) &  
     324                     & + e3u(ji,jj,jk,Krhs_a) * uu(ji,jj,jk,Krhs_a) * umask(ji,jj,jk) 
    371325               END DO 
    372326            END DO 
    373327            DO ji = 1, jpi 
    374                zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2) 
     328               zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 
    375329            END DO 
    376330                
    377331            DO jk = 1, jpkm1 
    378332               DO ji = 1, jpi 
    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) 
    381                END DO 
    382             END DO 
    383          ENDIF 
    384          ! 
    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  
     333                  uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) &  
     334                    & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 
     335               END DO 
     336            END DO 
     337         END DO 
    391338      ENDIF 
    392339      ! 
    … …  
    401348      !! 
    402349      INTEGER :: ji, jj 
     350      INTEGER :: istart, iend, jstart, jend 
    403351      !!----------------------------------------------------------------------   
    404352      ! 
    405353      IF( Agrif_Root() )   RETURN 
    406354      ! 
    407       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
     355      !--- West ---! 
     356      istart = 2 
     357      iend   = nbghostcells+1 
     358      DO ji = mi0(istart), mi1(iend) 
    408359         DO jj=1,jpj 
    409             va_e(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * hvr_e(2:nbghostcells+1,jj) 
    410             ! Specified fluxes: 
    411             ua_e(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * hur_e(2:nbghostcells+1,jj) 
    412             ! Characteristics method (only if ghostcells=1): 
    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)) ) 
    415          END DO 
    416       ENDIF 
    417       ! 
    418       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
     360            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     361            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     362         END DO 
     363      END DO 
     364      ! 
     365      !--- East ---! 
     366      istart = jpiglo-nbghostcells 
     367      iend   = jpiglo-1 
     368      DO ji = mi0(istart), mi1(iend) 
    419369         DO jj=1,jpj 
    420             va_e(nlci-nbghostcells:nlci-1,jj)   = vbdy_e(1:nbghostcells,jj) * hvr_e(nlci-nbghostcells:nlci-1,jj) 
    421             ! Specified fluxes: 
    422             ua_e(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * hur_e(nlci-nbghostcells-1:nlci-2,jj) 
    423             ! Characteristics method (only if ghostcells=1): 
    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)) ) 
    426          END DO 
    427       ENDIF 
    428       ! 
    429       IF((nbondj == -1).OR.(nbondj == 2)) THEN 
     370            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     371         END DO 
     372      END DO 
     373      istart = jpiglo-nbghostcells-1 
     374      iend   = jpiglo-2 
     375      DO ji = mi0(istart), mi1(iend) 
     376         DO jj=1,jpj 
     377            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     378         END DO 
     379      END DO 
     380      ! 
     381      !--- South ---! 
     382      jstart = 2 
     383      jend   = nbghostcells+1 
     384      DO jj = mj0(jstart), mj1(jend) 
    430385         DO ji=1,jpi 
    431             ua_e(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * hur_e(ji,2:nbghostcells+1) 
    432             ! Specified fluxes: 
    433             va_e(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * hvr_e(ji,2:nbghostcells+1) 
    434             ! Characteristics method (only if ghostcells=1): 
    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)) ) 
    437          END DO 
    438       ENDIF 
    439       ! 
    440       IF((nbondj == 1).OR.(nbondj == 2)) THEN 
     386            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     387            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     388         END DO 
     389      END DO 
     390      ! 
     391      !--- North ---! 
     392      jstart = jpjglo-nbghostcells 
     393      jend   = jpjglo-1 
     394      DO jj = mj0(jstart), mj1(jend) 
    441395         DO ji=1,jpi 
    442             ua_e(ji,nlcj-nbghostcells:nlcj-1)   = ubdy_n(ji,1:nbghostcells) * hur_e(ji,nlcj-nbghostcells:nlcj-1) 
    443             ! Specified fluxes: 
    444             va_e(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * hvr_e(ji,nlcj-nbghostcells-1:nlcj-2) 
    445             ! Characteristics method (only if ghostcells=1): 
    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)) ) 
    448          END DO 
    449       ENDIF 
     396            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     397         END DO 
     398      END DO 
     399      jstart = jpjglo-nbghostcells-1 
     400      jend   = jpjglo-2 
     401      DO jj = mj0(jstart), mj1(jend) 
     402         DO ji=1,jpi 
     403            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     404         END DO 
     405      END DO 
    450406      ! 
    451407   END SUBROUTINE Agrif_dyn_ts 
    452408 
     409   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv ) 
     410      !!---------------------------------------------------------------------- 
     411      !!                  ***  ROUTINE Agrif_dyn_ts_flux  *** 
     412      !!----------------------------------------------------------------------   
     413      INTEGER, INTENT(in) ::   jn 
     414      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv 
     415      !! 
     416      INTEGER :: ji, jj 
     417      INTEGER :: istart, iend, jstart, jend 
     418      !!----------------------------------------------------------------------   
     419      ! 
     420      IF( Agrif_Root() )   RETURN 
     421      ! 
     422      !--- West ---! 
     423      istart = 2 
     424      iend   = nbghostcells+1 
     425      DO ji = mi0(istart), mi1(iend) 
     426         DO jj=1,jpj 
     427            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     428            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     429         END DO 
     430      END DO 
     431      ! 
     432      !--- East ---! 
     433      istart = jpiglo-nbghostcells 
     434      iend   = jpiglo-1 
     435      DO ji = mi0(istart), mi1(iend) 
     436         DO jj=1,jpj 
     437            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     438         END DO 
     439      END DO 
     440      istart = jpiglo-nbghostcells-1 
     441      iend   = jpiglo-2 
     442      DO ji = mi0(istart), mi1(iend) 
     443         DO jj=1,jpj 
     444            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     445         END DO 
     446      END DO 
     447      ! 
     448      !--- South ---! 
     449      jstart = 2 
     450      jend   = nbghostcells+1 
     451      DO jj = mj0(jstart), mj1(jend) 
     452         DO ji=1,jpi 
     453            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     454            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     455         END DO 
     456      END DO 
     457      ! 
     458      !--- North ---! 
     459      jstart = jpjglo-nbghostcells 
     460      jend   = jpjglo-1 
     461      DO jj = mj0(jstart), mj1(jend) 
     462         DO ji=1,jpi 
     463            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     464         END DO 
     465      END DO 
     466      jstart = jpjglo-nbghostcells-1 
     467      jend   = jpjglo-2 
     468      DO jj = mj0(jstart), mj1(jend) 
     469         DO ji=1,jpi 
     470            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     471         END DO 
     472      END DO 
     473      ! 
     474   END SUBROUTINE Agrif_dyn_ts_flux 
    453475 
    454476   SUBROUTINE Agrif_dta_ts( kt ) 
    … …  
    470492      ! 
    471493      ! Interpolate barotropic fluxes 
    472       Agrif_SpecialValue=0._wp 
     494      Agrif_SpecialValue = 0._wp 
    473495      Agrif_UseSpecialValue = ln_spc_dyn 
     496      ! 
     497      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners) 
     498      utint_stage(:,:) = 0 
     499      vtint_stage(:,:) = 0 
    474500      ! 
    475501      IF( ll_int_cons ) THEN  ! Conservative interpolation 
    476502         ! order matters here !!!!!! 
    477503         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 ) 
     504         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
     505         ! 
    479506         bdy_tinterp = 1 
    480507         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After 
    481          CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
     508         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )   
     509         ! 
    482510         bdy_tinterp = 2 
    483511         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before 
    484          CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )          
     512         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )    
    485513      ELSE ! Linear interpolation 
    486          bdy_tinterp = 0 
    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 
     514         ! 
     515         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
    491516         CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 
    492517         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 
    … …  
    503528      INTEGER, INTENT(in) ::   kt 
    504529      ! 
    505       INTEGER  :: ji, jj, indx, indy 
     530      INTEGER  :: ji, jj 
     531      INTEGER  :: istart, iend, jstart, jend 
    506532      !!----------------------------------------------------------------------   
    507533      ! 
    … …  
    516542      ! 
    517543      ! --- West --- ! 
    518       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    519          indx = 1+nbghostcells 
     544      istart = 2 
     545      iend   = 1 + nbghostcells 
     546      DO ji = mi0(istart), mi1(iend) 
    520547         DO jj = 1, jpj 
    521             DO ji = 2, indx 
    522                ssha(ji,jj) = hbdy_w(ji-1,jj) 
    523             ENDDO 
     548            ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 
    524549         ENDDO 
    525       ENDIF 
     550      ENDDO 
    526551      ! 
    527552      ! --- East --- ! 
    528       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    529          indx = nlci-nbghostcells 
     553      istart = jpiglo - nbghostcells 
     554      iend   = jpiglo - 1 
     555      DO ji = mi0(istart), mi1(iend) 
    530556         DO jj = 1, jpj 
    531             DO ji = indx, nlci-1 
    532                ssha(ji,jj) = hbdy_e(ji-indx+1,jj) 
    533             ENDDO 
     557            ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 
    534558         ENDDO 
    535       ENDIF 
     559      ENDDO 
    536560      ! 
    537561      ! --- South --- ! 
    538       IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    539          indy = 1+nbghostcells 
    540          DO jj = 2, indy 
    541             DO ji = 1, jpi 
    542                ssha(ji,jj) = hbdy_s(ji,jj-1) 
    543             ENDDO 
     562      jstart = 2 
     563      jend   = 1 + nbghostcells 
     564      DO jj = mj0(jstart), mj1(jend) 
     565         DO ji = 1, jpi 
     566            ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 
    544567         ENDDO 
    545       ENDIF 
     568      ENDDO 
    546569      ! 
    547570      ! --- North --- ! 
    548       IF((nbondj == 1).OR.(nbondj == 2)) THEN 
    549          indy = nlcj-nbghostcells 
    550          DO jj = indy, nlcj-1 
    551             DO ji = 1, jpi 
    552                ssha(ji,jj) = hbdy_n(ji,jj-indy+1) 
    553             ENDDO 
     571      jstart = jpjglo - nbghostcells 
     572      jend   = jpjglo - 1 
     573      DO jj = mj0(jstart), mj1(jend) 
     574         DO ji = 1, jpi 
     575            ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 
    554576         ENDDO 
    555       ENDIF 
     577      ENDDO 
    556578      ! 
    557579   END SUBROUTINE Agrif_ssh 
    … …  
    564586      INTEGER, INTENT(in) ::   jn 
    565587      !! 
    566       INTEGER :: ji, jj, indx, indy 
    567       !!----------------------------------------------------------------------   
    568       !! clem ghost (starting at i,j=1 is important I think otherwise you introduce a grad(ssh)/=0 at point 2) 
     588      INTEGER :: ji, jj 
     589      INTEGER  :: istart, iend, jstart, jend 
     590      !!----------------------------------------------------------------------   
    569591      ! 
    570592      IF( Agrif_Root() )   RETURN 
    571593      ! 
    572594      ! --- West --- ! 
    573       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    574          indx = 1+nbghostcells 
     595      istart = 2 
     596      iend   = 1+nbghostcells 
     597      DO ji = mi0(istart), mi1(iend) 
    575598         DO jj = 1, jpj 
    576             DO ji = 2, indx 
    577                ssha_e(ji,jj) = hbdy_w(ji-1,jj) 
    578             ENDDO 
     599            ssha_e(ji,jj) = hbdy(ji,jj) 
    579600         ENDDO 
    580       ENDIF 
     601      ENDDO 
    581602      ! 
    582603      ! --- East --- ! 
    583       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    584          indx = nlci-nbghostcells 
     604      istart = jpiglo - nbghostcells 
     605      iend   = jpiglo - 1 
     606      DO ji = mi0(istart), mi1(iend) 
    585607         DO jj = 1, jpj 
    586             DO ji = indx, nlci-1 
    587                ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj) 
    588             ENDDO 
     608            ssha_e(ji,jj) = hbdy(ji,jj) 
    589609         ENDDO 
    590       ENDIF 
     610      ENDDO 
    591611      ! 
    592612      ! --- South --- ! 
    593       IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    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 
     613      jstart = 2 
     614      jend   = 1+nbghostcells 
     615      DO jj = mj0(jstart), mj1(jend) 
     616         DO ji = 1, jpi 
     617            ssha_e(ji,jj) = hbdy(ji,jj) 
    599618         ENDDO 
    600       ENDIF 
     619      ENDDO 
    601620      ! 
    602621      ! --- North --- ! 
    603       IF((nbondj == 1).OR.(nbondj == 2)) THEN 
    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 
     622      jstart = jpjglo - nbghostcells 
     623      jend   = jpjglo - 1 
     624      DO jj = mj0(jstart), mj1(jend) 
     625         DO ji = 1, jpi 
     626            ssha_e(ji,jj) = hbdy(ji,jj) 
    609627         ENDDO 
    610       ENDIF 
     628      ENDDO 
    611629      ! 
    612630   END SUBROUTINE Agrif_ssh_ts 
    … …  
    634652    
    635653 
    636    SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     654   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    637655      !!---------------------------------------------------------------------- 
    638656      !!                  *** ROUTINE interptsn *** 
    … …  
    641659      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    642660      LOGICAL                                     , INTENT(in   ) ::   before 
    643       INTEGER                                     , INTENT(in   ) ::   nb , ndir 
    644       ! 
    645       INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices 
    646       INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
    647       REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7 
    648       LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     661      ! 
     662      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices 
     663      INTEGER  ::   N_in, N_out 
    649664      ! 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 
     665      REAL(wp) :: zhtot 
     666      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 
    652667      REAL(wp), DIMENSION(k1:k2) :: h_in 
    653668      REAL(wp), DIMENSION(1:jpk) :: h_out 
    654       REAL(wp) :: h_diff 
     669      !!---------------------------------------------------------------------- 
    655670 
    656671      IF( before ) THEN          
    … …  
    659674               DO jj=j1,j2 
    660675                 DO ji=i1,i2 
    661                        ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     676                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) 
    662677                 END DO 
    663678              END DO 
    … …  
    666681 
    667682# if defined key_vertical 
     683        ! Interpolate thicknesses 
     684        ! Warning: these are masked, hence extrapolated prior interpolation. 
    668685        DO jk=k1,k2 
    669686           DO jj=j1,j2 
    670687              DO ji=i1,i2 
    671                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     688                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    672689              END DO 
    673690           END DO 
    674691        END DO 
     692 
     693        ! Extrapolate thicknesses in partial bottom cells: 
     694        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     695        IF (ln_zps) THEN 
     696           DO jj=j1,j2 
     697              DO ji=i1,i2 
     698                  jk = mbkt(ji,jj) 
     699                  ptab(ji,jj,jk,jpts+1) = 0._wp 
     700              END DO 
     701           END DO            
     702        END IF 
     703      
     704        ! Save ssh at last level: 
     705        IF (.NOT.ln_linssh) THEN 
     706           ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
     707        ELSE 
     708           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
     709        END IF       
    675710# endif 
    676711      ELSE  
    677712 
    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) 
    680  
    681 # if defined key_vertical               
     713# if defined key_vertical  
     714         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp  
     715             
    682716         DO jj=j1,j2 
    683717            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 
     718               ts(ji,jj,:,:,Krhs_a) = 0._wp 
     719               N_in = mbkt_parent(ji,jj) 
     720               zhtot = 0._wp 
     721               DO jk=1,N_in !k2 = jpk of parent grid 
     722                  IF (jk==N_in) THEN 
     723                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 
     724                  ELSE 
     725                     h_in(jk) = ptab(ji,jj,jk,n2) 
     726                  ENDIF 
     727                  zhtot = zhtot + h_in(jk) 
    694728                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    695                   h_in(N_in) = ptab(ji,jj,jk,n2) 
    696729               END DO 
    697730               N_out = 0 
    698731               DO jk=1,jpk ! jpk of child grid 
    699                   IF (tmask(iref,jref,jk) == 0) EXIT  
     732                  IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    700733                  N_out = N_out + 1 
    701                   h_out(jk) = e3t_n(iref,jref,jk) 
     734                  h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
    702735               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 
     736               IF (N_in*N_out > 0) THEN 
     737                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts) 
    707738               ENDIF 
    708739            ENDDO 
    709740         ENDDO 
    710741# else 
    711          ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts) 
    712 # endif 
    713742         ! 
    714743         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 
    717  
    718          IF ( .NOT.lk_agrif_clp ) THEN  
    719             ! 
    720             imin = i1 ; imax = i2 
    721             jmin = j1 ; jmax = j2 
    722             !  
    723             ! Remove CORNERS 
    724             IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells 
    725             IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1 
    726             IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells 
    727             IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1       
    728             ! 
    729             IF( eastern_side ) THEN 
    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                ! 
    737                ibdy = nlci-nbghostcells 
    738                DO jn = 1, jpts 
    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) 
    740                   DO jk = 1, jpkm1 
    741                      DO jj = jmin,jmax 
    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) 
    744                         ELSE 
    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) 
    749                            ENDIF 
    750                         ENDIF 
    751                      END DO 
    752                   END DO 
    753                   ! Restore ghost points: 
    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) 
    755                END DO 
    756             ENDIF 
    757             !  
    758             IF( northern_side ) THEN 
    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                ! 
    766                jbdy = nlcj-nbghostcells          
    767                DO jn = 1, jpts 
    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) 
    769                   DO jk = 1, jpkm1 
    770                      DO ji = imin,imax 
    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) 
    773                         ELSE 
    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) 
    778                            ENDIF 
    779                         ENDIF 
    780                      END DO 
    781                   END DO 
    782                   ! Restore ghost points: 
    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) 
    784                END DO 
    785             ENDIF 
    786             ! 
    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                !     
    795                ibdy = 1+nbghostcells        
    796                DO jn = 1, jpts 
    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) 
    798                   DO jk = 1, jpkm1 
    799                      DO jj = jmin,jmax 
    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) 
    802                         ELSE 
    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 
    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) 
    807                            ENDIF 
    808                         ENDIF 
    809                      END DO 
    810                   END DO 
    811                   ! Restore ghost points: 
    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) 
    813                END DO 
    814             ENDIF 
    815             ! 
    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                !   
    824                jbdy=1+nbghostcells         
    825                DO jn = 1, jpts 
    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 
    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) 
    831                         ELSE 
    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 
    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) 
    836                            ENDIF 
    837                         ENDIF 
    838                      END DO 
    839                   END DO 
    840                   ! Restore ghost points: 
    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) 
    842                END DO 
    843             ENDIF 
    844             ! 
    845          ENDIF 
     744            ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
     745         END DO 
     746# endif 
     747 
    846748      ENDIF 
    847749      ! 
    848750   END SUBROUTINE interptsn 
    849751 
    850    SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     752   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before ) 
    851753      !!---------------------------------------------------------------------- 
    852754      !!                  ***  ROUTINE interpsshn  *** 
    … …  
    855757      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    856758      LOGICAL                         , INTENT(in   ) ::   before 
    857       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    858       ! 
    859       LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     759      ! 
    860760      !!----------------------------------------------------------------------   
    861761      ! 
    862762      IF( before) THEN 
    863          ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2) 
     763         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a) 
    864764      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) 
    869          !! clem ghost 
    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) 
     765         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
    874766      ENDIF 
    875767      ! 
    876768   END SUBROUTINE interpsshn 
    877769 
    878    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     770   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    879771      !!---------------------------------------------------------------------- 
    880772      !!                  *** ROUTINE interpun *** 
    … …  
    884776      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    885777      LOGICAL, INTENT(in) :: before 
    886       INTEGER, INTENT(in) :: nb , ndir 
    887778      !! 
    888779      INTEGER :: ji,jj,jk 
    889       REAL(wp) :: zrhoy 
     780      REAL(wp) :: zrhoy, zhtot 
    890781      ! vertical interpolation: 
    891782      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    892783      REAL(wp), DIMENSION(1:jpk) :: h_out 
    893       INTEGER  :: N_in, N_out, iref 
     784      INTEGER  :: N_in, N_out 
    894785      REAL(wp) :: h_diff 
    895       LOGICAL  :: western_side, eastern_side 
    896786      !!---------------------------------------------     
    897787      ! 
    … …  
    900790            DO jj=j1,j2 
    901791               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))  
     792                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk))  
    903793# if defined key_vertical 
    904                   ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 
     794                  ! Interpolate thicknesses (masked for subsequent extrapolation) 
     795                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    905796# endif 
    906797               END DO 
    907798            END DO 
    908799         END DO 
     800# if defined key_vertical 
     801         ! Extrapolate thicknesses in partial bottom cells: 
     802         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     803         IF (ln_zps) THEN 
     804            DO jj=j1,j2 
     805               DO ji=i1,i2 
     806                  jk = mbku(ji,jj) 
     807                  ptab(ji,jj,jk,2) = 0._wp 
     808               END DO 
     809            END DO            
     810         END IF 
     811        ! Save ssh at last level: 
     812        ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     813        IF (.NOT.ln_linssh) THEN 
     814           ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
     815           DO jk=1,jpk 
     816              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 
     817           END DO 
     818           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 
     819        END IF  
     820# endif 
     821         ! 
    909822      ELSE 
    910823         zrhoy = Agrif_rhoy() 
    911824# if defined key_vertical 
    912825! VERTICAL REFINEMENT BEGIN 
    913          western_side  = (nb == 1).AND.(ndir == 1) 
    914          eastern_side  = (nb == 1).AND.(ndir == 2) 
     826 
     827         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
    915828 
    916829         DO ji=i1,i2 
    917             iref = ji 
    918             IF (western_side) iref = MAX(2,ji) 
    919             IF (eastern_side) iref = MIN(nlci-2,ji) 
    920830            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)  
     831               uu(ji,jj,:,Krhs_a) = 0._wp 
     832               N_in = mbku_parent(ji,jj) 
     833               zhtot = 0._wp 
     834               DO jk=1,N_in 
     835                  IF (jk==N_in) THEN 
     836                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     837                  ELSE 
     838                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     839                  ENDIF 
     840                  zhtot = zhtot + h_in(jk) 
     841                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 
    927842              ENDDO 
    928           
    929               IF (N_in == 0) THEN 
    930                  ua(ji,jj,:) = 0._wp 
    931                  CYCLE 
    932               ENDIF 
    933           
     843                   
    934844              N_out = 0 
    935845              DO jk=1,jpk 
    936                  if (umask(iref,jj,jk) == 0) EXIT 
     846                 if (umask(ji,jj,jk) == 0) EXIT 
    937847                 N_out = N_out + 1 
    938                  h_out(N_out) = e3u_a(iref,jj,jk) 
     848                 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 
    939849              ENDDO 
    940           
    941               IF (N_out == 0) THEN 
    942                  ua(ji,jj,:) = 0._wp 
    943                  CYCLE 
     850              IF (N_in*N_out > 0) THEN 
     851                 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 
    944852              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) 
    955853            ENDDO 
    956854         ENDDO 
    … …  
    959857         DO jk = 1, jpkm1 
    960858            DO jj=j1,j2 
    961                ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 
     859               uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) ) 
    962860            END DO 
    963861         END DO 
    … …  
    968866   END SUBROUTINE interpun 
    969867 
    970    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     868   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    971869      !!---------------------------------------------------------------------- 
    972870      !!                  *** ROUTINE interpvn *** 
    … …  
    976874      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    977875      LOGICAL, INTENT(in) :: before 
    978       INTEGER, INTENT(in) :: nb , ndir 
    979876      ! 
    980877      INTEGER :: ji,jj,jk 
    … …  
    983880      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    984881      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 
     882      INTEGER  :: N_in, N_out 
     883      REAL(wp) :: h_diff, zhtot 
    988884      !!---------------------------------------------     
    989885      !       
    … …  
    992888            DO jj=j1,j2 
    993889               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)) 
     890                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk)) 
    995891# if defined key_vertical 
    996                   ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 
     892                  ! Interpolate thicknesses (masked for subsequent extrapolation) 
     893                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    997894# endif 
    998895               END DO 
    999896            END DO 
    1000897         END DO 
     898# if defined key_vertical 
     899         ! Extrapolate thicknesses in partial bottom cells: 
     900         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     901         IF (ln_zps) THEN 
     902            DO jj=j1,j2 
     903               DO ji=i1,i2 
     904                  jk = mbkv(ji,jj) 
     905                  ptab(ji,jj,jk,2) = 0._wp 
     906               END DO 
     907            END DO            
     908         END IF 
     909        ! Save ssh at last level: 
     910        ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     911        IF (.NOT.ln_linssh) THEN 
     912           ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
     913           DO jk=1,jpk 
     914              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 
     915           END DO 
     916           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 
     917        END IF  
     918# endif 
    1001919      ELSE        
    1002920         zrhox = Agrif_rhox() 
    1003921# if defined key_vertical 
    1004922 
    1005          southern_side = (nb == 2).AND.(ndir == 1) 
    1006          northern_side = (nb == 2).AND.(ndir == 2) 
     923         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
    1007924 
    1008925         DO jj=j1,j2 
    1009             jref = jj 
    1010             IF (southern_side) jref = MAX(2,jj) 
    1011             IF (northern_side) jref = MIN(nlcj-2,jj) 
    1012926            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 
     927               vv(ji,jj,:,Krhs_a) = 0._wp 
     928               N_in = mbkv_parent(ji,jj) 
     929               zhtot = 0._wp 
     930               DO jk=1,N_in 
     931                  IF (jk==N_in) THEN 
     932                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     933                  ELSE 
     934                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     935                  ENDIF 
     936                  zhtot = zhtot + h_in(jk) 
     937                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 
     938              ENDDO 
    1024939          
    1025940               N_out = 0 
    1026941               DO jk=1,jpk 
    1027                   if (vmask(ji,jref,jk) == 0) EXIT 
     942                  if (vmask(ji,jj,jk) == 0) EXIT 
    1028943                  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 
     944                  h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 
     945               END DO 
     946               IF (N_in*N_out > 0) THEN 
     947                  call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 
    1034948               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) 
    1036949            END DO 
    1037950         END DO 
    1038951# else 
    1039952         DO jk = 1, jpkm1 
    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) ) 
     953            vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) ) 
    1041954         END DO 
    1042955# endif 
    … …  
    1045958   END SUBROUTINE interpvn 
    1046959 
    1047    SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     960   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before) 
    1048961      !!---------------------------------------------------------------------- 
    1049962      !!                  ***  ROUTINE interpunb  *** 
    … …  
    1052965      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    1053966      LOGICAL                         , INTENT(in   ) ::   before 
    1054       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    1055967      ! 
    1056968      INTEGER  ::   ji, jj 
    1057969      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff 
    1058       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    1059970      !!----------------------------------------------------------------------   
    1060971      ! 
    1061972      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) 
     973         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * uu_b(i1:i2,j1:j2,Kmm_a) 
    1063974      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) 
    1068975         zrhoy = Agrif_Rhoy() 
    1069976         zrhot = Agrif_rhot() 
    … …  
    1071978         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot 
    1072979         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)        & 
    1076                &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
    1077          ELSEIF( bdy_tinterp == 2 ) THEN 
    1078             ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
    1079                &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
    1080          ELSE 
    1081             ztcoeff = 1 
    1082          ENDIF 
    1083          !    
    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)  
    1088          !             
    1089          IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 
    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) 
    1094          ENDIF 
    1095       ENDIF 
     980         !  
     981         DO ji = i1, i2 
     982            DO jj = j1, j2 
     983               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     984                  IF    ( utint_stage(ji,jj) == 1  ) THEN 
     985                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     986                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     987                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN 
     988                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     989                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     990                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN                 
     991                     ztcoeff = 1._wp 
     992                  ELSE 
     993                     ztcoeff = 0._wp 
     994                  ENDIF 
     995                  !    
     996                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     997                  !             
     998                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN 
     999                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1) 
     1000                  ENDIF 
     1001                  ! 
     1002                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1 
     1003               ENDIF 
     1004            END DO 
     1005         END DO 
     1006      END IF 
    10961007      !  
    10971008   END SUBROUTINE interpunb 
    10981009 
    10991010 
    1100    SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     1011   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before ) 
    11011012      !!---------------------------------------------------------------------- 
    11021013      !!                  ***  ROUTINE interpvnb  *** 
    … …  
    11051016      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    11061017      LOGICAL                         , INTENT(in   ) ::   before 
    1107       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    1108       ! 
    1109       INTEGER  ::   ji,jj 
     1018      ! 
     1019      INTEGER  ::   ji, jj 
    11101020      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff    
    1111       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    11121021      !!----------------------------------------------------------------------   
    11131022      !  
    11141023      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) 
     1024         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a) 
    11161025      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) 
    11211026         zrhox = Agrif_Rhox() 
    11221027         zrhot = Agrif_rhot() 
    11231028         ! Time indexes bounds for integration 
    11241029         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)        & 
    1128                &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
    1129          ELSEIF( bdy_tinterp == 2 ) THEN 
    1130             ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
    1131                &               - zt0        * (       zt0 - 1._wp)**2._wp )  
    1132          ELSE 
    1133             ztcoeff = 1 
    1134          ENDIF 
    1135          !! clem ghost 
    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)  
    1140          !             
    1141          IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN 
    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) 
    1146          ENDIF 
     1030         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot  
     1031         !      
     1032         DO ji = i1, i2 
     1033            DO jj = j1, j2 
     1034               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     1035                  IF    ( vtint_stage(ji,jj) == 1  ) THEN 
     1036                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     1037                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     1038                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN 
     1039                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     1040                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     1041                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                 
     1042                     ztcoeff = 1._wp 
     1043                  ELSE 
     1044                     ztcoeff = 0._wp 
     1045                  ENDIF 
     1046                  !    
     1047                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     1048                  !             
     1049                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN 
     1050                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1) 
     1051                  ENDIF 
     1052                  ! 
     1053                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1 
     1054               ENDIF 
     1055            END DO 
     1056         END DO           
    11471057      ENDIF 
    11481058      ! 
    … …  
    11501060 
    11511061 
    1152    SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     1062   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before ) 
    11531063      !!---------------------------------------------------------------------- 
    11541064      !!                  ***  ROUTINE interpub2b  *** 
    … …  
    11571067      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    11581068      LOGICAL                         , INTENT(in   ) ::   before 
    1159       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    11601069      ! 
    11611070      INTEGER  ::   ji,jj 
    1162       REAL(wp) ::   zrhot, zt0, zt1,zat 
    1163       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
     1071      REAL(wp) ::   zrhot, zt0, zt1, zat 
    11641072      !!----------------------------------------------------------------------   
    11651073      IF( before ) THEN 
    … …  
    11701078         ENDIF 
    11711079      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: 
    1181          zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    & 
    1182             &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    )  
    1183          !! clem ghost 
    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)  
    1188       ENDIF 
    1189       !  
    1190    END SUBROUTINE interpub2b 
    1191     
    1192  
    1193    SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir ) 
    1194       !!---------------------------------------------------------------------- 
    1195       !!                  ***  ROUTINE interpvb2b  *** 
    1196       !!----------------------------------------------------------------------   
    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 
    1205       !!----------------------------------------------------------------------   
    1206       ! 
    1207       IF( before ) THEN 
    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 
    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) 
    12181080         zrhot = Agrif_rhot() 
    12191081         ! Time indexes bounds for integration 
    … …  
    12241086            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    )  
    12251087         ! 
    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)  
     1088         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)  
     1089         ! 
     1090         ! Update interpolation stage: 
     1091         utint_stage(i1:i2,j1:j2) = 1 
     1092      ENDIF 
     1093      !  
     1094   END SUBROUTINE interpub2b 
     1095    
     1096 
     1097   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before ) 
     1098      !!---------------------------------------------------------------------- 
     1099      !!                  ***  ROUTINE interpvb2b  *** 
     1100      !!----------------------------------------------------------------------   
     1101      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1102      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1103      LOGICAL                         , INTENT(in   ) ::   before 
     1104      ! 
     1105      INTEGER ::   ji,jj 
     1106      REAL(wp) ::   zrhot, zt0, zt1, zat 
     1107      !!----------------------------------------------------------------------   
     1108      ! 
     1109      IF( before ) THEN 
     1110         IF ( ln_bt_fw ) THEN 
     1111            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1112         ELSE 
     1113            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1114         ENDIF 
     1115      ELSE       
     1116         zrhot = Agrif_rhot() 
     1117         ! Time indexes bounds for integration 
     1118         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot 
     1119         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
     1120         ! Polynomial interpolation coefficients: 
     1121         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    & 
     1122            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    )  
     1123         ! 
     1124         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
     1125         ! 
     1126         ! update interpolation stage: 
     1127         vtint_stage(i1:i2,j1:j2) = 1 
    12301128      ENDIF 
    12311129      !       
    … …  
    12331131 
    12341132 
    1235    SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 
     1133   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 
    12361134      !!---------------------------------------------------------------------- 
    12371135      !!                  ***  ROUTINE interpe3t  *** 
    … …  
    12401138      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
    12411139      LOGICAL                              , INTENT(in   ) :: before 
    1242       INTEGER                              , INTENT(in   ) :: nb , ndir 
    12431140      ! 
    12441141      INTEGER :: ji, jj, jk 
    1245       LOGICAL :: western_side, eastern_side, northern_side, southern_side 
    12461142      !!----------------------------------------------------------------------   
    12471143      !     
    … …  
    12491145         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 
    12501146      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) 
    12551147         ! 
    12561148         DO jk = k1, k2 
    12571149            DO jj = j1, j2 
    12581150               DO ji = i1, i2 
    1259                   ! 
    12601151                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN 
    1261                      IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN 
    1262                         WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 
    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 
    1266                         WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk 
    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 
    1270                         WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 
    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 
    1274                         WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk 
    1275                         WRITE(numout,*)  ptab(ji,jj,jk), e3t_0(ji,jj,jk) 
    1276                         kindic_agr = kindic_agr + 1 
    1277                      ENDIF 
     1152                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  &  
     1153                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), & 
     1154                     &                 ji+nimpp-1, jj+njmpp-1, jk 
     1155                     kindic_agr = kindic_agr + 1 
    12781156                  ENDIF 
    12791157               END DO 
    … …  
    12841162      !  
    12851163   END SUBROUTINE interpe3t 
    1286  
    1287  
    1288    SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 
    1289       !!---------------------------------------------------------------------- 
    1290       !!                  ***  ROUTINE interpumsk  *** 
    1291       !!----------------------------------------------------------------------   
    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 
    1296       ! 
    1297       INTEGER ::   ji, jj, jk 
    1298       LOGICAL ::   western_side, eastern_side    
    1299       !!----------------------------------------------------------------------   
    1300       !     
    1301       IF( before ) THEN 
    1302          ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2) 
    1303       ELSE 
    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 
    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 
    1323          END DO 
    1324          ! 
    1325       ENDIF 
    1326       !  
    1327    END SUBROUTINE interpumsk 
    1328  
    1329  
    1330    SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 
    1331       !!---------------------------------------------------------------------- 
    1332       !!                  ***  ROUTINE interpvmsk  *** 
    1333       !!----------------------------------------------------------------------   
    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 
    1338       ! 
    1339       INTEGER ::   ji, jj, jk 
    1340       LOGICAL ::   northern_side, southern_side      
    1341       !!----------------------------------------------------------------------   
    1342       !     
    1343       IF( before ) THEN 
    1344          ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2) 
    1345       ELSE 
    1346          southern_side = (nb == 2).AND.(ndir == 1) 
    1347          northern_side = (nb == 2).AND.(ndir == 2) 
    1348          DO jk = k1, k2 
    1349             DO jj = j1, j2 
    1350                DO ji = i1, i2 
    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 
    1365          END DO 
    1366          ! 
    1367       ENDIF 
    1368       !  
    1369    END SUBROUTINE interpvmsk 
    13701164 
    13711165 
    … …  
    13771171      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab 
    13781172      LOGICAL                                    , INTENT(in   ) ::   before 
    1379       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    1380       REAL(wp), DIMENSION(1:jpk) :: h_out 
    1381       INTEGER  :: N_in, N_out, ji, jj, jk 
     1173      ! 
     1174      INTEGER  :: ji, jj, jk 
     1175      INTEGER  :: N_in, N_out 
     1176      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in 
     1177      REAL(wp), DIMENSION(1:jpk) :: z_out 
    13821178      !!----------------------------------------------------------------------   
    13831179      !       
    … …  
    13901186           END DO 
    13911187        END DO 
    1392 #ifdef key_vertical          
     1188 
     1189# if defined key_vertical 
     1190        ! Interpolate thicknesses 
     1191        ! Warning: these are masked, hence extrapolated prior interpolation. 
    13931192        DO jk=k1,k2 
    13941193           DO jj=j1,j2 
    13951194              DO ji=i1,i2 
    1396                  ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk)  
     1195                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    13971196              END DO 
    13981197           END DO 
    13991198        END DO 
    1400 #endif 
     1199 
     1200        ! Extrapolate thicknesses in partial bottom cells: 
     1201        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     1202        IF (ln_zps) THEN 
     1203           DO jj=j1,j2 
     1204              DO ji=i1,i2 
     1205                  jk = mbkt(ji,jj) 
     1206                  ptab(ji,jj,jk,2) = 0._wp 
     1207              END DO 
     1208           END DO            
     1209        END IF 
     1210      
     1211        ! Save ssh at last level: 
     1212        IF (.NOT.ln_linssh) THEN 
     1213           ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
     1214        ELSE 
     1215           ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     1216        END IF       
     1217# endif 
    14011218      ELSE  
    14021219#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) 
    1411                   h_in(N_in) = ptab(ji,jj,jk,2) 
    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) 
     1220         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
     1221         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 
     1222             
     1223         DO jj = j1, j2 
     1224            DO ji =i1, i2 
     1225               N_in = mbkt_parent(ji,jj) 
     1226               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 
     1227               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 
     1228               DO jk = N_in, 1, -1  ! Parent vertical grid                
     1229                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 
     1230                    tabin(jk) = ptab(ji,jj,jk,1) 
     1231               END DO 
     1232               N_out = mbkt(ji,jj)  
     1233               DO jk = 1, N_out        ! Child vertical grid 
     1234                  z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 
    14181235               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) 
     1236               IF (N_in*N_out > 0) THEN 
     1237                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 
    14211238               ENDIF 
    14221239            ENDDO 
    … …  
    14281245      ! 
    14291246   END SUBROUTINE interpavm 
     1247 
     1248# if defined key_vertical 
     1249   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 
     1250      !!---------------------------------------------------------------------- 
     1251      !!                  ***  ROUTINE interpsshn  *** 
     1252      !!----------------------------------------------------------------------   
     1253      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1254      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1255      LOGICAL                         , INTENT(in   ) ::   before 
     1256      ! 
     1257      !!----------------------------------------------------------------------   
     1258      ! 
     1259      IF( before) THEN 
     1260         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp) 
     1261      ELSE 
     1262         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2)) 
     1263      ENDIF 
     1264      ! 
     1265   END SUBROUTINE interpmbkt 
     1266 
     1267   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before ) 
     1268      !!---------------------------------------------------------------------- 
     1269      !!                  ***  ROUTINE interpsshn  *** 
     1270      !!----------------------------------------------------------------------   
     1271      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1272      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1273      LOGICAL                         , INTENT(in   ) ::   before 
     1274      ! 
     1275      !!----------------------------------------------------------------------   
     1276      ! 
     1277      IF( before) THEN 
     1278         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) 
     1279      ELSE 
     1280         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     1281      ENDIF 
     1282      ! 
     1283   END SUBROUTINE interpht0 
     1284#endif 
    14301285 
    14311286#else 
Note: See TracChangeset for help on using the changeset viewer.