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 12123 for NEMO/branches/2019/dev_AGRIF-01-05_merged/src/NST/agrif_oce_interp.F90 – NEMO

Ignore:
Timestamp:
2019-12-09T13:55:34+01:00 (4 years ago)
Author:
jchanut
Message:

Merge devs from #2199 and #2222 in trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_AGRIF-01-05_merged/src/NST/agrif_oce_interp.F90

    r10068 r12123  
    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 
     
    7881      ! 
    7982      INTEGER ::   ji, jj, jk       ! dummy loop indices 
    80       INTEGER ::   j1, j2, i1, i2 
    8183      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2 
    8284      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb 
     
    9395      Agrif_UseSpecialValue = .FALSE. 
    9496      ! 
    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  
    10397      ! --- 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 
     98      ibdy1 = 2 
     99      ibdy2 = 1+nbghostcells  
     100      ! 
     101      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     102         DO ji = mi0(ibdy1), mi1(ibdy2) 
     103            ua_b(ji,:) = 0._wp 
     104 
    110105            DO jk = 1, jpkm1 
    111106               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 
     107                  ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
     108               END DO 
     109            END DO 
     110 
    116111            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 
     112               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj) 
     113            END DO 
     114         END DO 
     115      ENDIF 
     116      ! 
     117      DO ji = mi0(ibdy1), mi1(ibdy2) 
     118         zub(ji,:) = 0._wp    ! Correct transport 
    130119         DO jk = 1, jpkm1 
    131120            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) 
     121               zub(ji,jj) = zub(ji,jj) &  
     122                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk)*umask(ji,jj,jk) 
    134123            END DO 
    135124         END DO 
    136125         DO jj=1,jpj 
    137             zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
     126            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 
    138127         END DO 
    139128             
    140129         DO jk = 1, jpkm1 
    141130            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 
     131               ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk) 
     132            END DO 
     133         END DO 
     134      END DO 
    146135             
    147          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    148             zvb(ibdy1:ibdy2,:) = 0._wp 
     136      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     137         DO ji = mi0(ibdy1), mi1(ibdy2) 
     138            zvb(ji,:) = 0._wp 
    149139            DO jk = 1, jpkm1 
    150140               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) 
     141                  zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    153142               END DO 
    154143            END DO 
    155144            DO jj = 1, jpj 
    156                zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj) 
     145               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 
    157146            END DO 
    158147            DO jk = 1, jpkm1 
    159148               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  
     149                  va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk) 
     150               END DO 
     151            END DO 
     152         END DO 
    172153      ENDIF 
    173154 
    174155      ! --- 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 
     156      ibdy1 = jpiglo-1-nbghostcells 
     157      ibdy2 = jpiglo-2  
     158      ! 
     159      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     160         DO ji = mi0(ibdy1), mi1(ibdy2) 
     161            ua_b(ji,:) = 0._wp 
    181162            DO jk = 1, jpkm1 
    182163               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) 
     164                  ua_b(ji,jj) = ua_b(ji,jj) &  
     165                      & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    185166               END DO 
    186167            END DO 
    187168            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 
     169               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj) 
     170            END DO 
     171         END DO 
     172      ENDIF 
     173      ! 
     174      DO ji = mi0(ibdy1), mi1(ibdy2) 
     175         zub(ji,:) = 0._wp    ! Correct transport 
    201176         DO jk = 1, jpkm1 
    202177            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) 
     178               zub(ji,jj) = zub(ji,jj) &  
     179                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk) * umask(ji,jj,jk) 
    205180            END DO 
    206181         END DO 
    207182         DO jj=1,jpj 
    208             zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj) 
     183            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 
    209184         END DO 
    210185             
    211186         DO jk = 1, jpkm1 
    212187            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 
     188               ua(ji,jj,jk) = ( ua(ji,jj,jk) &  
     189                 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk) 
     190            END DO 
     191         END DO 
     192      END DO 
    217193             
    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 
     194      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     195         ibdy1 = jpiglo-nbghostcells 
     196         ibdy2 = jpiglo-1  
     197         DO ji = mi0(ibdy1), mi1(ibdy2) 
     198            zvb(ji,:) = 0._wp 
    222199            DO jk = 1, jpkm1 
    223200               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) 
     201                  zvb(ji,jj) = zvb(ji,jj) & 
     202                     & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    226203               END DO 
    227204            END DO 
    228205            DO jj = 1, jpj 
    229                zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj) 
     206               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 
    230207            END DO 
    231208            DO jk = 1, jpkm1 
    232209               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  
     210                  va(ji,jj,jk) = ( va(ji,jj,jk) &  
     211                      & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk) 
     212               END DO 
     213            END DO 
     214         END DO 
    245215      ENDIF 
    246216 
    247217      ! --- 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 
     218      jbdy1 = 2 
     219      jbdy2 = 1+nbghostcells  
     220      ! 
     221      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     222         DO jj = mj0(jbdy1), mj1(jbdy2) 
     223            va_b(:,jj) = 0._wp 
    254224            DO jk = 1, jpkm1 
    255225               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) 
     226                  va_b(ji,jj) = va_b(ji,jj) &  
     227                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    258228               END DO 
    259229            END DO 
    260230            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 
     231               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)      
     232            END DO 
     233         END DO 
     234      ENDIF 
     235      ! 
     236      DO jj = mj0(jbdy1), mj1(jbdy2) 
     237         zvb(:,jj) = 0._wp    ! Correct transport 
    274238         DO jk=1,jpkm1 
    275239            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) 
     240               zvb(ji,jj) = zvb(ji,jj) &  
     241                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    278242            END DO 
    279243         END DO 
    280244         DO ji = 1, jpi 
    281             zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
     245            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 
    282246         END DO 
    283247 
    284248         DO jk = 1, jpkm1 
    285249            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 
     250               va(ji,jj,jk) = ( va(ji,jj,jk) &  
     251                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk) 
     252            END DO 
     253         END DO 
     254      END DO 
    290255             
    291          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
    292             zub(:,jbdy1:jbdy2) = 0._wp 
     256      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     257         DO jj = mj0(jbdy1), mj1(jbdy2) 
     258            zub(:,jj) = 0._wp 
    293259            DO jk = 1, jpkm1 
    294260               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) 
     261                  zub(ji,jj) = zub(ji,jj) &  
     262                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    297263               END DO 
    298264            END DO 
    299265            DO ji = 1, jpi 
    300                zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2) 
     266               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 
    301267            END DO 
    302268                
    303269            DO jk = 1, jpkm1 
    304270               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  
     271                  ua(ji,jj,jk) = ( ua(ji,jj,jk) &  
     272                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk) 
     273               END DO 
     274            END DO 
     275         END DO 
    317276      ENDIF 
    318277 
    319278      ! --- 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 
     279      jbdy1 = jpjglo-1-nbghostcells 
     280      jbdy2 = jpjglo-2  
     281      ! 
     282      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport 
     283         DO jj = mj0(jbdy1), mj1(jbdy2) 
     284            va_b(:,jj) = 0._wp 
    326285            DO jk = 1, jpkm1 
    327286               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) 
     287                  va_b(ji,jj) = va_b(ji,jj) &  
     288                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    330289               END DO 
    331290            END DO 
    332291            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 
     292               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj) 
     293            END DO 
     294         END DO 
     295      ENDIF 
     296      ! 
     297      DO jj = mj0(jbdy1), mj1(jbdy2) 
     298         zvb(:,jj) = 0._wp    ! Correct transport 
    346299         DO jk=1,jpkm1 
    347300            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) 
     301               zvb(ji,jj) = zvb(ji,jj) &  
     302                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 
    350303            END DO 
    351304         END DO 
    352305         DO ji = 1, jpi 
    353             zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2) 
     306            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj) 
    354307         END DO 
    355308 
    356309         DO jk = 1, jpkm1 
    357310            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 
     311               va(ji,jj,jk) = ( va(ji,jj,jk) &  
     312                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk) 
     313            END DO 
     314         END DO 
     315      END DO 
    362316             
    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 
     317      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate 
     318         jbdy1 = jpjglo-nbghostcells 
     319         jbdy2 = jpjglo-1 
     320         DO jj = mj0(jbdy1), mj1(jbdy2) 
     321            zub(:,jj) = 0._wp 
    367322            DO jk = 1, jpkm1 
    368323               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) 
     324                  zub(ji,jj) = zub(ji,jj) &  
     325                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 
    371326               END DO 
    372327            END DO 
    373328            DO ji = 1, jpi 
    374                zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2) 
     329               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj) 
    375330            END DO 
    376331                
    377332            DO jk = 1, jpkm1 
    378333               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  
     334                  ua(ji,jj,jk) = ( ua(ji,jj,jk) &  
     335                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk) 
     336               END DO 
     337            END DO 
     338         END DO 
    391339      ENDIF 
    392340      ! 
     
    401349      !! 
    402350      INTEGER :: ji, jj 
     351      INTEGER :: istart, iend, jstart, jend 
    403352      !!----------------------------------------------------------------------   
    404353      ! 
    405354      IF( Agrif_Root() )   RETURN 
    406355      ! 
    407       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
     356      !--- West ---! 
     357      istart = 2 
     358      iend   = nbghostcells+1 
     359      DO ji = mi0(istart), mi1(iend) 
    408360         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 
     361            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     362            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     363         END DO 
     364      END DO 
     365      ! 
     366      !--- East ---! 
     367      istart = jpiglo-nbghostcells 
     368      iend   = jpiglo-1 
     369      DO ji = mi0(istart), mi1(iend) 
    419370         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 
     371            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     372         END DO 
     373      END DO 
     374      istart = jpiglo-nbghostcells-1 
     375      iend   = jpiglo-2 
     376      DO ji = mi0(istart), mi1(iend) 
     377         DO jj=1,jpj 
     378            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     379         END DO 
     380      END DO 
     381      ! 
     382      !--- South ---! 
     383      jstart = 2 
     384      jend   = nbghostcells+1 
     385      DO jj = mj0(jstart), mj1(jend) 
    430386         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 
     387            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     388            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     389         END DO 
     390      END DO 
     391      ! 
     392      !--- North ---! 
     393      jstart = jpjglo-nbghostcells 
     394      jend   = jpjglo-1 
     395      DO jj = mj0(jstart), mj1(jend) 
    441396         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 
     397            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) 
     398         END DO 
     399      END DO 
     400      jstart = jpjglo-nbghostcells-1 
     401      jend   = jpjglo-2 
     402      DO jj = mj0(jstart), mj1(jend) 
     403         DO ji=1,jpi 
     404            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) 
     405         END DO 
     406      END DO 
    450407      ! 
    451408   END SUBROUTINE Agrif_dyn_ts 
    452409 
     410   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv ) 
     411      !!---------------------------------------------------------------------- 
     412      !!                  ***  ROUTINE Agrif_dyn_ts_flux  *** 
     413      !!----------------------------------------------------------------------   
     414      INTEGER, INTENT(in) ::   jn 
     415      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv 
     416      !! 
     417      INTEGER :: ji, jj 
     418      INTEGER :: istart, iend, jstart, jend 
     419      !!----------------------------------------------------------------------   
     420      ! 
     421      IF( Agrif_Root() )   RETURN 
     422      ! 
     423      !--- West ---! 
     424      istart = 2 
     425      iend   = nbghostcells+1 
     426      DO ji = mi0(istart), mi1(iend) 
     427         DO jj=1,jpj 
     428            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     429            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     430         END DO 
     431      END DO 
     432      ! 
     433      !--- East ---! 
     434      istart = jpiglo-nbghostcells 
     435      iend   = jpiglo-1 
     436      DO ji = mi0(istart), mi1(iend) 
     437         DO jj=1,jpj 
     438            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     439         END DO 
     440      END DO 
     441      istart = jpiglo-nbghostcells-1 
     442      iend   = jpiglo-2 
     443      DO ji = mi0(istart), mi1(iend) 
     444         DO jj=1,jpj 
     445            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     446         END DO 
     447      END DO 
     448      ! 
     449      !--- South ---! 
     450      jstart = 2 
     451      jend   = nbghostcells+1 
     452      DO jj = mj0(jstart), mj1(jend) 
     453         DO ji=1,jpi 
     454            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     455            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     456         END DO 
     457      END DO 
     458      ! 
     459      !--- North ---! 
     460      jstart = jpjglo-nbghostcells 
     461      jend   = jpjglo-1 
     462      DO jj = mj0(jstart), mj1(jend) 
     463         DO ji=1,jpi 
     464            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj) 
     465         END DO 
     466      END DO 
     467      jstart = jpjglo-nbghostcells-1 
     468      jend   = jpjglo-2 
     469      DO jj = mj0(jstart), mj1(jend) 
     470         DO ji=1,jpi 
     471            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj) 
     472         END DO 
     473      END DO 
     474      ! 
     475   END SUBROUTINE Agrif_dyn_ts_flux 
    453476 
    454477   SUBROUTINE Agrif_dta_ts( kt ) 
     
    470493      ! 
    471494      ! Interpolate barotropic fluxes 
    472       Agrif_SpecialValue=0._wp 
     495      Agrif_SpecialValue = 0._wp 
    473496      Agrif_UseSpecialValue = ln_spc_dyn 
     497      ! 
     498      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners) 
     499      utint_stage(:,:) = 0 
     500      vtint_stage(:,:) = 0 
    474501      ! 
    475502      IF( ll_int_cons ) THEN  ! Conservative interpolation 
    476503         ! order matters here !!!!!! 
    477504         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 ) 
     505         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
     506         ! 
    479507         bdy_tinterp = 1 
    480508         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After 
    481          CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
     509         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )   
     510         ! 
    482511         bdy_tinterp = 2 
    483512         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before 
    484          CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )          
     513         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )    
    485514      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 
     515         ! 
     516         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp  
    491517         CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 
    492518         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 
     
    503529      INTEGER, INTENT(in) ::   kt 
    504530      ! 
    505       INTEGER  :: ji, jj, indx, indy 
     531      INTEGER  :: ji, jj 
     532      INTEGER  :: istart, iend, jstart, jend 
    506533      !!----------------------------------------------------------------------   
    507534      ! 
     
    516543      ! 
    517544      ! --- West --- ! 
    518       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    519          indx = 1+nbghostcells 
     545      istart = 2 
     546      iend   = 1 + nbghostcells 
     547      DO ji = mi0(istart), mi1(iend) 
    520548         DO jj = 1, jpj 
    521             DO ji = 2, indx 
    522                ssha(ji,jj) = hbdy_w(ji-1,jj) 
    523             ENDDO 
     549            ssha(ji,jj) = hbdy(ji,jj) 
    524550         ENDDO 
    525       ENDIF 
     551      ENDDO 
    526552      ! 
    527553      ! --- East --- ! 
    528       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    529          indx = nlci-nbghostcells 
     554      istart = jpiglo - nbghostcells 
     555      iend   = jpiglo - 1 
     556      DO ji = mi0(istart), mi1(iend) 
    530557         DO jj = 1, jpj 
    531             DO ji = indx, nlci-1 
    532                ssha(ji,jj) = hbdy_e(ji-indx+1,jj) 
    533             ENDDO 
     558            ssha(ji,jj) = hbdy(ji,jj) 
    534559         ENDDO 
    535       ENDIF 
     560      ENDDO 
    536561      ! 
    537562      ! --- 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 
     563      jstart = 2 
     564      jend   = 1 + nbghostcells 
     565      DO jj = mj0(jstart), mj1(jend) 
     566         DO ji = 1, jpi 
     567            ssha(ji,jj) = hbdy(ji,jj) 
    544568         ENDDO 
    545       ENDIF 
     569      ENDDO 
    546570      ! 
    547571      ! --- 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 
     572      jstart = jpjglo - nbghostcells 
     573      jend   = jpjglo - 1 
     574      DO jj = mj0(jstart), mj1(jend) 
     575         DO ji = 1, jpi 
     576            ssha(ji,jj) = hbdy(ji,jj) 
    554577         ENDDO 
    555       ENDIF 
     578      ENDDO 
    556579      ! 
    557580   END SUBROUTINE Agrif_ssh 
     
    564587      INTEGER, INTENT(in) ::   jn 
    565588      !! 
    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) 
     589      INTEGER :: ji, jj 
     590      INTEGER  :: istart, iend, jstart, jend 
     591      !!----------------------------------------------------------------------   
    569592      ! 
    570593      IF( Agrif_Root() )   RETURN 
    571594      ! 
    572595      ! --- West --- ! 
    573       IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    574          indx = 1+nbghostcells 
     596      istart = 2 
     597      iend   = 1+nbghostcells 
     598      DO ji = mi0(istart), mi1(iend) 
    575599         DO jj = 1, jpj 
    576             DO ji = 2, indx 
    577                ssha_e(ji,jj) = hbdy_w(ji-1,jj) 
    578             ENDDO 
     600            ssha_e(ji,jj) = hbdy(ji,jj) 
    579601         ENDDO 
    580       ENDIF 
     602      ENDDO 
    581603      ! 
    582604      ! --- East --- ! 
    583       IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    584          indx = nlci-nbghostcells 
     605      istart = jpiglo - nbghostcells 
     606      iend   = jpiglo - 1 
     607      DO ji = mi0(istart), mi1(iend) 
    585608         DO jj = 1, jpj 
    586             DO ji = indx, nlci-1 
    587                ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj) 
    588             ENDDO 
     609            ssha_e(ji,jj) = hbdy(ji,jj) 
    589610         ENDDO 
    590       ENDIF 
     611      ENDDO 
    591612      ! 
    592613      ! --- 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 
     614      jstart = 2 
     615      jend   = 1+nbghostcells 
     616      DO jj = mj0(jstart), mj1(jend) 
     617         DO ji = 1, jpi 
     618            ssha_e(ji,jj) = hbdy(ji,jj) 
    599619         ENDDO 
    600       ENDIF 
     620      ENDDO 
    601621      ! 
    602622      ! --- 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 
     623      jstart = jpjglo - nbghostcells 
     624      jend   = jpjglo - 1 
     625      DO jj = mj0(jstart), mj1(jend) 
     626         DO ji = 1, jpi 
     627            ssha_e(ji,jj) = hbdy(ji,jj) 
    609628         ENDDO 
    610       ENDIF 
     629      ENDDO 
    611630      ! 
    612631   END SUBROUTINE Agrif_ssh_ts 
     
    634653    
    635654 
    636    SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     655   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    637656      !!---------------------------------------------------------------------- 
    638657      !!                  *** ROUTINE interptsn *** 
     
    641660      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    642661      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 
     662      ! 
     663      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices 
     664      INTEGER  ::   N_in, N_out 
    649665      ! 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 
     666      REAL(wp) :: zhtot 
     667      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 
    652668      REAL(wp), DIMENSION(k1:k2) :: h_in 
    653669      REAL(wp), DIMENSION(1:jpk) :: h_out 
    654       REAL(wp) :: h_diff 
     670      !!---------------------------------------------------------------------- 
    655671 
    656672      IF( before ) THEN          
     
    666682 
    667683# if defined key_vertical 
     684        ! Interpolate thicknesses 
     685        ! Warning: these are masked, hence extrapolated prior interpolation. 
    668686        DO jk=k1,k2 
    669687           DO jj=j1,j2 
    670688              DO ji=i1,i2 
    671                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     689                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
    672690              END DO 
    673691           END DO 
    674692        END DO 
     693 
     694        ! Extrapolate thicknesses in partial bottom cells: 
     695        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     696        IF (ln_zps) THEN 
     697           DO jj=j1,j2 
     698              DO ji=i1,i2 
     699                  jk = mbkt(ji,jj) 
     700                  ptab(ji,jj,jk,jpts+1) = 0._wp 
     701              END DO 
     702           END DO            
     703        END IF 
     704      
     705        ! Save ssh at last level: 
     706        IF (.NOT.ln_linssh) THEN 
     707           ptab(i1:i2,j1:j2,k2,jpts+1) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)  
     708        ELSE 
     709           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
     710        END IF       
    675711# endif 
    676712      ELSE  
    677713 
    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               
     714# if defined key_vertical  
     715         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp  
     716             
    682717         DO jj=j1,j2 
    683718            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 
     719               tsa(ji,jj,:,:) = 0._wp 
     720               N_in = mbkt_parent(ji,jj) 
     721               zhtot = 0._wp 
     722               DO jk=1,N_in !k2 = jpk of parent grid 
     723                  IF (jk==N_in) THEN 
     724                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 
     725                  ELSE 
     726                     h_in(jk) = ptab(ji,jj,jk,n2) 
     727                  ENDIF 
     728                  zhtot = zhtot + h_in(jk) 
    694729                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    695                   h_in(N_in) = ptab(ji,jj,jk,n2) 
    696730               END DO 
    697731               N_out = 0 
    698732               DO jk=1,jpk ! jpk of child grid 
    699                   IF (tmask(iref,jref,jk) == 0) EXIT  
     733                  IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    700734                  N_out = N_out + 1 
    701                   h_out(jk) = e3t_n(iref,jref,jk) 
     735                  h_out(jk) = e3t_a(ji,jj,jk) 
    702736               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 
     737               IF (N_in*N_out > 0) THEN 
     738                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tsa(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    707739               ENDIF 
    708740            ENDDO 
    709741         ENDDO 
    710742# else 
    711          ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts) 
    712 # endif 
    713743         ! 
    714744         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 
     745            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
     746         END DO 
     747# endif 
     748 
    846749      ENDIF 
    847750      ! 
    848751   END SUBROUTINE interptsn 
    849752 
    850    SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     753   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before ) 
    851754      !!---------------------------------------------------------------------- 
    852755      !!                  ***  ROUTINE interpsshn  *** 
     
    855758      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    856759      LOGICAL                         , INTENT(in   ) ::   before 
    857       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    858       ! 
    859       LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     760      ! 
    860761      !!----------------------------------------------------------------------   
    861762      ! 
     
    863764         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2) 
    864765      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) 
     766         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
    874767      ENDIF 
    875768      ! 
    876769   END SUBROUTINE interpsshn 
    877770 
    878    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     771   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    879772      !!---------------------------------------------------------------------- 
    880773      !!                  *** ROUTINE interpun *** 
     
    884777      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    885778      LOGICAL, INTENT(in) :: before 
    886       INTEGER, INTENT(in) :: nb , ndir 
    887779      !! 
    888780      INTEGER :: ji,jj,jk 
    889       REAL(wp) :: zrhoy 
     781      REAL(wp) :: zrhoy, zhtot 
    890782      ! vertical interpolation: 
    891783      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    892784      REAL(wp), DIMENSION(1:jpk) :: h_out 
    893       INTEGER  :: N_in, N_out, iref 
     785      INTEGER  :: N_in, N_out 
    894786      REAL(wp) :: h_diff 
    895       LOGICAL  :: western_side, eastern_side 
    896787      !!---------------------------------------------     
    897788      ! 
     
    902793                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk))  
    903794# if defined key_vertical 
    904                   ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 
     795                  ! Interpolate thicknesses (masked for subsequent extrapolation) 
     796                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) 
    905797# endif 
    906798               END DO 
    907799            END DO 
    908800         END DO 
     801# if defined key_vertical 
     802         ! Extrapolate thicknesses in partial bottom cells: 
     803         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     804         IF (ln_zps) THEN 
     805            DO jj=j1,j2 
     806               DO ji=i1,i2 
     807                  jk = mbku(ji,jj) 
     808                  ptab(ji,jj,jk,2) = 0._wp 
     809               END DO 
     810            END DO            
     811         END IF 
     812        ! Save ssh at last level: 
     813        ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     814        IF (.NOT.ln_linssh) THEN 
     815           ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
     816           DO jk=1,jpk 
     817              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u_n(i1:i2,j1:j2,jk) * umask(i1:i2,j1:j2,jk) 
     818           END DO 
     819           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 
     820        END IF  
     821# endif 
     822         ! 
    909823      ELSE 
    910824         zrhoy = Agrif_rhoy() 
    911825# if defined key_vertical 
    912826! VERTICAL REFINEMENT BEGIN 
    913          western_side  = (nb == 1).AND.(ndir == 1) 
    914          eastern_side  = (nb == 1).AND.(ndir == 2) 
     827 
     828         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
    915829 
    916830         DO ji=i1,i2 
    917             iref = ji 
    918             IF (western_side) iref = MAX(2,ji) 
    919             IF (eastern_side) iref = MIN(nlci-2,ji) 
    920831            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)  
     832               ua(ji,jj,:) = 0._wp 
     833               N_in = mbku_parent(ji,jj) 
     834               zhtot = 0._wp 
     835               DO jk=1,N_in 
     836                  IF (jk==N_in) THEN 
     837                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     838                  ELSE 
     839                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     840                  ENDIF 
     841                  zhtot = zhtot + h_in(jk) 
     842                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 
    927843              ENDDO 
    928           
    929               IF (N_in == 0) THEN 
    930                  ua(ji,jj,:) = 0._wp 
    931                  CYCLE 
    932               ENDIF 
    933           
     844                   
    934845              N_out = 0 
    935846              DO jk=1,jpk 
    936                  if (umask(iref,jj,jk) == 0) EXIT 
     847                 if (umask(ji,jj,jk) == 0) EXIT 
    937848                 N_out = N_out + 1 
    938                  h_out(N_out) = e3u_a(iref,jj,jk) 
     849                 h_out(N_out) = e3u_a(ji,jj,jk) 
    939850              ENDDO 
    940           
    941               IF (N_out == 0) THEN 
    942                  ua(ji,jj,:) = 0._wp 
    943                  CYCLE 
     851              IF (N_in*N_out > 0) THEN 
     852                 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,1) 
    944853              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) 
    955854            ENDDO 
    956855         ENDDO 
     
    968867   END SUBROUTINE interpun 
    969868 
    970    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     869   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    971870      !!---------------------------------------------------------------------- 
    972871      !!                  *** ROUTINE interpvn *** 
     
    976875      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    977876      LOGICAL, INTENT(in) :: before 
    978       INTEGER, INTENT(in) :: nb , ndir 
    979877      ! 
    980878      INTEGER :: ji,jj,jk 
     
    983881      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    984882      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 
     883      INTEGER  :: N_in, N_out 
     884      REAL(wp) :: h_diff, zhtot 
    988885      !!---------------------------------------------     
    989886      !       
     
    994891                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 
    995892# if defined key_vertical 
     893                  ! Interpolate thicknesses (masked for subsequent extrapolation) 
    996894                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 
    997895# endif 
     
    999897            END DO 
    1000898         END DO 
     899# if defined key_vertical 
     900         ! Extrapolate thicknesses in partial bottom cells: 
     901         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     902         IF (ln_zps) THEN 
     903            DO jj=j1,j2 
     904               DO ji=i1,i2 
     905                  jk = mbkv(ji,jj) 
     906                  ptab(ji,jj,jk,2) = 0._wp 
     907               END DO 
     908            END DO            
     909         END IF 
     910        ! Save ssh at last level: 
     911        ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     912        IF (.NOT.ln_linssh) THEN 
     913           ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
     914           DO jk=1,jpk 
     915              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v_n(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk) 
     916           END DO 
     917           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 
     918        END IF  
     919# endif 
    1001920      ELSE        
    1002921         zrhox = Agrif_rhox() 
    1003922# if defined key_vertical 
    1004923 
    1005          southern_side = (nb == 2).AND.(ndir == 1) 
    1006          northern_side = (nb == 2).AND.(ndir == 2) 
     924         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
    1007925 
    1008926         DO jj=j1,j2 
    1009             jref = jj 
    1010             IF (southern_side) jref = MAX(2,jj) 
    1011             IF (northern_side) jref = MIN(nlcj-2,jj) 
    1012927            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 
     928               va(ji,jj,:) = 0._wp 
     929               N_in = mbkv_parent(ji,jj) 
     930               zhtot = 0._wp 
     931               DO jk=1,N_in 
     932                  IF (jk==N_in) THEN 
     933                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 
     934                  ELSE 
     935                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)  
     936                  ENDIF 
     937                  zhtot = zhtot + h_in(jk) 
     938                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 
     939              ENDDO 
    1024940          
    1025941               N_out = 0 
    1026942               DO jk=1,jpk 
    1027                   if (vmask(ji,jref,jk) == 0) EXIT 
     943                  if (vmask(ji,jj,jk) == 0) EXIT 
    1028944                  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 
     945                  h_out(N_out) = e3v_a(ji,jj,jk) 
     946               END DO 
     947               IF (N_in*N_out > 0) THEN 
     948                  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,1) 
    1034949               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) 
    1036950            END DO 
    1037951         END DO 
     
    1045959   END SUBROUTINE interpvn 
    1046960 
    1047    SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     961   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before) 
    1048962      !!---------------------------------------------------------------------- 
    1049963      !!                  ***  ROUTINE interpunb  *** 
     
    1052966      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    1053967      LOGICAL                         , INTENT(in   ) ::   before 
    1054       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    1055968      ! 
    1056969      INTEGER  ::   ji, jj 
    1057970      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff 
    1058       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    1059971      !!----------------------------------------------------------------------   
    1060972      ! 
     
    1062974         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2) 
    1063975      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) 
    1068976         zrhoy = Agrif_Rhoy() 
    1069977         zrhot = Agrif_rhot() 
     
    1071979         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot 
    1072980         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 
     981         !  
     982         DO ji = i1, i2 
     983            DO jj = j1, j2 
     984               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     985                  IF    ( utint_stage(ji,jj) == 1  ) THEN 
     986                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     987                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     988                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN 
     989                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     990                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     991                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN                 
     992                     ztcoeff = 1._wp 
     993                  ELSE 
     994                     ztcoeff = 0._wp 
     995                  ENDIF 
     996                  !    
     997                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     998                  !             
     999                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN 
     1000                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1) 
     1001                  ENDIF 
     1002                  ! 
     1003                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1 
     1004               ENDIF 
     1005            END DO 
     1006         END DO 
     1007      END IF 
    10961008      !  
    10971009   END SUBROUTINE interpunb 
    10981010 
    10991011 
    1100    SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     1012   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before ) 
    11011013      !!---------------------------------------------------------------------- 
    11021014      !!                  ***  ROUTINE interpvnb  *** 
     
    11051017      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    11061018      LOGICAL                         , INTENT(in   ) ::   before 
    1107       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    1108       ! 
    1109       INTEGER  ::   ji,jj 
     1019      ! 
     1020      INTEGER  ::   ji, jj 
    11101021      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff    
    1111       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    11121022      !!----------------------------------------------------------------------   
    11131023      !  
     
    11151025         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2) 
    11161026      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) 
    11211027         zrhox = Agrif_Rhox() 
    11221028         zrhot = Agrif_rhot() 
    11231029         ! Time indexes bounds for integration 
    11241030         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 
     1031         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot  
     1032         !      
     1033         DO ji = i1, i2 
     1034            DO jj = j1, j2 
     1035               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     1036                  IF    ( vtint_stage(ji,jj) == 1  ) THEN 
     1037                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     1038                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     1039                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN 
     1040                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     1041                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     1042                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                 
     1043                     ztcoeff = 1._wp 
     1044                  ELSE 
     1045                     ztcoeff = 0._wp 
     1046                  ENDIF 
     1047                  !    
     1048                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     1049                  !             
     1050                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN 
     1051                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1) 
     1052                  ENDIF 
     1053                  ! 
     1054                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1 
     1055               ENDIF 
     1056            END DO 
     1057         END DO           
    11471058      ENDIF 
    11481059      ! 
     
    11501061 
    11511062 
    1152    SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir ) 
     1063   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before ) 
    11531064      !!---------------------------------------------------------------------- 
    11541065      !!                  ***  ROUTINE interpub2b  *** 
     
    11571068      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    11581069      LOGICAL                         , INTENT(in   ) ::   before 
    1159       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    11601070      ! 
    11611071      INTEGER  ::   ji,jj 
    1162       REAL(wp) ::   zrhot, zt0, zt1,zat 
    1163       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
     1072      REAL(wp) ::   zrhot, zt0, zt1, zat 
    11641073      !!----------------------------------------------------------------------   
    11651074      IF( before ) THEN 
     
    11701079         ENDIF 
    11711080      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) 
    12181081         zrhot = Agrif_rhot() 
    12191082         ! Time indexes bounds for integration 
     
    12241087            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    )  
    12251088         ! 
    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)  
     1089         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)  
     1090         ! 
     1091         ! Update interpolation stage: 
     1092         utint_stage(i1:i2,j1:j2) = 1 
     1093      ENDIF 
     1094      !  
     1095   END SUBROUTINE interpub2b 
     1096    
     1097 
     1098   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before ) 
     1099      !!---------------------------------------------------------------------- 
     1100      !!                  ***  ROUTINE interpvb2b  *** 
     1101      !!----------------------------------------------------------------------   
     1102      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1103      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1104      LOGICAL                         , INTENT(in   ) ::   before 
     1105      ! 
     1106      INTEGER ::   ji,jj 
     1107      REAL(wp) ::   zrhot, zt0, zt1, zat 
     1108      !!----------------------------------------------------------------------   
     1109      ! 
     1110      IF( before ) THEN 
     1111         IF ( ln_bt_fw ) THEN 
     1112            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1113         ELSE 
     1114            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1115         ENDIF 
     1116      ELSE       
     1117         zrhot = Agrif_rhot() 
     1118         ! Time indexes bounds for integration 
     1119         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot 
     1120         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
     1121         ! Polynomial interpolation coefficients: 
     1122         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    & 
     1123            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    )  
     1124         ! 
     1125         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
     1126         ! 
     1127         ! update interpolation stage: 
     1128         vtint_stage(i1:i2,j1:j2) = 1 
    12301129      ENDIF 
    12311130      !       
     
    12331132 
    12341133 
    1235    SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 
     1134   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before ) 
    12361135      !!---------------------------------------------------------------------- 
    12371136      !!                  ***  ROUTINE interpe3t  *** 
     
    12401139      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
    12411140      LOGICAL                              , INTENT(in   ) :: before 
    1242       INTEGER                              , INTENT(in   ) :: nb , ndir 
    12431141      ! 
    12441142      INTEGER :: ji, jj, jk 
    1245       LOGICAL :: western_side, eastern_side, northern_side, southern_side 
    12461143      !!----------------------------------------------------------------------   
    12471144      !     
     
    12491146         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 
    12501147      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) 
    12551148         ! 
    12561149         DO jk = k1, k2 
    12571150            DO jj = j1, j2 
    12581151               DO ji = i1, i2 
    1259                   ! 
    12601152                  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 
     1153                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  &  
     1154                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), & 
     1155                     &                 ji+nimpp-1, jj+njmpp-1, jk 
     1156                     kindic_agr = kindic_agr + 1 
    12781157                  ENDIF 
    12791158               END DO 
     
    12841163      !  
    12851164   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 
    13701165 
    13711166 
     
    13771172      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab 
    13781173      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 
     1174      ! 
     1175      INTEGER  :: ji, jj, jk 
     1176      INTEGER  :: N_in, N_out 
     1177      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in 
     1178      REAL(wp), DIMENSION(1:jpk) :: z_out 
    13821179      !!----------------------------------------------------------------------   
    13831180      !       
     
    13901187           END DO 
    13911188        END DO 
    1392 #ifdef key_vertical          
     1189 
     1190# if defined key_vertical 
     1191        ! Interpolate thicknesses 
     1192        ! Warning: these are masked, hence extrapolated prior interpolation. 
    13931193        DO jk=k1,k2 
    13941194           DO jj=j1,j2 
    13951195              DO ji=i1,i2 
    1396                  ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk)  
     1196                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
    13971197              END DO 
    13981198           END DO 
    13991199        END DO 
    1400 #endif 
     1200 
     1201        ! Extrapolate thicknesses in partial bottom cells: 
     1202        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     1203        IF (ln_zps) THEN 
     1204           DO jj=j1,j2 
     1205              DO ji=i1,i2 
     1206                  jk = mbkt(ji,jj) 
     1207                  ptab(ji,jj,jk,2) = 0._wp 
     1208              END DO 
     1209           END DO            
     1210        END IF 
     1211      
     1212        ! Save ssh at last level: 
     1213        IF (.NOT.ln_linssh) THEN 
     1214           ptab(i1:i2,j1:j2,k2,2) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)  
     1215        ELSE 
     1216           ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     1217        END IF       
     1218# endif 
    14011219      ELSE  
    14021220#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) 
     1221         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
     1222         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 
     1223             
     1224         DO jj = j1, j2 
     1225            DO ji =i1, i2 
     1226               N_in = mbkt_parent(ji,jj) 
     1227               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 
     1228               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 
     1229               DO jk = N_in, 1, -1  ! Parent vertical grid                
     1230                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 
     1231                    tabin(jk) = ptab(ji,jj,jk,1) 
    14121232               END DO 
    14131233               N_out = 0 
    1414                DO jk=1,jpk ! jpk of child grid 
    1415                   IF (wmask(ji,jj,jk) == 0) EXIT  
     1234          z_out(1) = 0._wp  
     1235               DO jk = 2, jpk       ! Child vertical grid 
     1236                  IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    14161237                  N_out = N_out + 1 
    1417                   h_out(jk) = e3t_n(ji,jj,jk) 
     1238                  z_out(jk) = z_out(jk-1) + e3t_n(ji,jj,jk-1) 
    14181239               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) 
     1240               IF (N_in*N_out > 0) THEN 
     1241                  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) 
    14211242               ENDIF 
    14221243            ENDDO 
     
    14281249      ! 
    14291250   END SUBROUTINE interpavm 
     1251 
     1252# if defined key_vertical 
     1253   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 
     1254      !!---------------------------------------------------------------------- 
     1255      !!                  ***  ROUTINE interpsshn  *** 
     1256      !!----------------------------------------------------------------------   
     1257      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1258      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1259      LOGICAL                         , INTENT(in   ) ::   before 
     1260      ! 
     1261      !!----------------------------------------------------------------------   
     1262      ! 
     1263      IF( before) THEN 
     1264         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp) 
     1265      ELSE 
     1266         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2)) 
     1267      ENDIF 
     1268      ! 
     1269   END SUBROUTINE interpmbkt 
     1270 
     1271   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before ) 
     1272      !!---------------------------------------------------------------------- 
     1273      !!                  ***  ROUTINE interpsshn  *** 
     1274      !!----------------------------------------------------------------------   
     1275      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1276      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1277      LOGICAL                         , INTENT(in   ) ::   before 
     1278      ! 
     1279      !!----------------------------------------------------------------------   
     1280      ! 
     1281      IF( before) THEN 
     1282         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) 
     1283      ELSE 
     1284         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     1285      ENDIF 
     1286      ! 
     1287   END SUBROUTINE interpht0 
     1288#endif 
    14301289 
    14311290#else 
Note: See TracChangeset for help on using the changeset viewer.