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 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/NST/agrif_oce_sponge.F90 – NEMO

Ignore:
Timestamp:
2021-05-05T13:18:04+02:00 (3 years ago)
Author:
mcastril
Message:

[2021/HPC-11_mcastril_HPDAonline_DiagGPU] Update externals

Location:
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
         5^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8^/vendors/PPR@HEAD            ext/PPR 
        89 
        910# SETTE 
        10 ^/utils/CI/sette@13559        sette 
         11^/utils/CI/sette@14244        sette 

  • NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/NST/agrif_oce_sponge.F90

    r13312 r14789  
    44   !!====================================================================== 
    55   !!                   ***  MODULE  agrif_oce_interp  *** 
    6    !! AGRIF: sponge package for the ocean dynamics (OPA) 
     6   !! AGRIF: sponge package for the ocean dynamics (OCE) 
    77   !!====================================================================== 
    88   !! History :  2.0  !  2002-06  (XXX)  Original cade 
    … …  
    2828   PRIVATE 
    2929 
    30    PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
     30   PUBLIC Agrif_Sponge, Agrif_Sponge_2d, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
    3131   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge 
     32   PUBLIC interpunb_sponge, interpvnb_sponge 
    3233 
    3334   !! * Substitutions 
     35#  include "domzgr_substitute.h90" 
    3436#  include "do_loop_substitute.h90" 
    3537   !!---------------------------------------------------------------------- 
    … …  
    4547      !!---------------------------------------------------------------------- 
    4648      REAL(wp) ::   zcoef   ! local scalar 
    47        
     49      INTEGER  :: istart, iend, jstart, jend  
    4850      !!---------------------------------------------------------------------- 
    4951      ! 
    … …  
    5254      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
    5355 
    54       CALL Agrif_Sponge 
    5556      Agrif_SpecialValue    = 0._wp 
    5657      Agrif_UseSpecialValue = .TRUE. 
     58      l_vremap              = ln_vert_remap 
    5759      tabspongedone_tsn     = .FALSE. 
    5860      ! 
    59       CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge ) 
     61      CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge ) 
    6062      ! 
    6163      Agrif_UseSpecialValue = .FALSE. 
     64      l_vremap              = .FALSE. 
    6265#endif 
    6366      ! 
    … …  
    8083      Agrif_SpecialValue    = 0._wp 
    8184      Agrif_UseSpecialValue = ln_spc_dyn 
     85      l_vremap              = ln_vert_remap 
    8286      use_sign_north        = .TRUE. 
    8387      sign_north            = -1._wp 
    … …  
    9094      tabspongedone_v = .FALSE. 
    9195      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge ) 
     96       
     97      IF ( nn_shift_bar>0 ) THEN ! then split sponge between 2d and 3d 
     98         zcoef = REAL(Agrif_NbStepint(),wp)/REAL(Agrif_rhot()) ! forward tsplit 
     99         tabspongedone_u = .FALSE. 
     100         tabspongedone_v = .FALSE.          
     101         CALL Agrif_Bc_Variable( unb_sponge_id, calledweight=zcoef, procname=interpunb_sponge ) 
     102         ! 
     103         tabspongedone_u = .FALSE. 
     104         tabspongedone_v = .FALSE. 
     105         CALL Agrif_Bc_Variable( vnb_sponge_id, calledweight=zcoef, procname=interpvnb_sponge ) 
     106      ENDIF 
    92107      ! 
    93108      Agrif_UseSpecialValue = .FALSE. 
    94109      use_sign_north        = .FALSE. 
     110      l_vremap              = .FALSE. 
     111      ! 
    95112#endif 
    96113      ! 
    … …  
    109126      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
    110127      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
    111 #if defined key_vertical 
    112       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu 
    113       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv 
    114 #endif 
    115       REAL(wp), DIMENSION(jpjmax)  :: zmskwest,  zmskeast 
    116       REAL(wp), DIMENSION(jpimax)  :: zmsknorth, zmsksouth 
    117128      !!---------------------------------------------------------------------- 
    118129      ! 
    … …  
    123134      !|            |           |           |           | 
    124135      !fine :     t u t u t u t u t u t u t u t u t u t u t 
    125       !sponge val:0   0   0   1  5/6 4/6 3/6 2/6 1/6  0   0 
     136      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0  
    126137      !           |   ghost     | <-- sponge area  -- > | 
    127138      !           |   points    |                       | 
    … …  
    129140 
    130141#if defined SPONGE || defined SPONGE_TOP 
    131       IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
    132          ! 
    133          ! Retrieve masks at open boundaries: 
    134  
    135          IF( lk_west ) THEN                             ! --- West --- ! 
    136             ztabramp(:,:) = 0._wp 
    137             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    138             DO ji = mi0(ind1), mi1(ind1)                 
    139                ztabramp(ji,:) = ssumask(ji,:) 
    140             END DO 
    141             zmskwest(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1) 
    142             zmskwest(jpj+1:jpjmax) = 0._wp 
    143          ENDIF 
    144          IF( lk_east ) THEN                             ! --- East --- ! 
    145             ztabramp(:,:) = 0._wp 
    146             ind1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    147             DO ji = mi0(ind1), mi1(ind1)                  
    148                ztabramp(ji,:) = ssumask(ji,:) 
    149             END DO 
    150             zmskeast(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1) 
    151             zmskeast(jpj+1:jpjmax) = 0._wp 
    152          ENDIF 
    153          IF( lk_south ) THEN                            ! --- South --- ! 
    154             ztabramp(:,:) = 0._wp 
    155             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    156             DO jj = mj0(ind1), mj1(ind1)                  
    157                ztabramp(:,jj) = ssvmask(:,jj) 
    158             END DO 
    159             zmsksouth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2) 
    160             zmsksouth(jpi+1:jpimax) = 0._wp 
    161          ENDIF 
    162          IF( lk_north ) THEN                            ! --- North --- ! 
    163             ztabramp(:,:) = 0._wp 
    164             ind1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells 
    165             DO jj = mj0(ind1), mj1(ind1)                  
    166                ztabramp(:,jj) = ssvmask(:,jj) 
    167             END DO 
    168             zmsknorth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2) 
    169             zmsknorth(jpi+1:jpimax) = 0._wp 
    170          ENDIF 
    171  
    172          ! JC: SPONGE MASKING TO BE SORTED OUT: 
    173          zmskwest(:)  = 1._wp 
    174          zmskeast(:)  = 1._wp 
    175          zmsksouth(:) = 1._wp 
    176          zmsknorth(:) = 1._wp 
    177 #if defined key_mpp_mpi 
    178 !         CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 
    179 !         CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 
    180 !         CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 
    181 !         CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 
     142      ! Define ramp from boundaries towards domain interior at F-points 
     143      ! Store it in ztabramp 
     144 
     145      ispongearea    = nn_sponge_len * Agrif_irhox() 
     146      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     147      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     148      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     149          
     150      ztabramp(:,:) = 0._wp 
     151 
     152      IF( lk_west ) THEN                             ! --- West --- ! 
     153         ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     154         ind2 = nn_hls + 1 + nbghostcells + ispongearea  
     155         DO ji = mi0(ind1), mi1(ind2)    
     156            DO jj = 1, jpj                
     157               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     158            END DO 
     159         END DO 
     160         ! ghost cells: 
     161         ind1 = 1 
     162         ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
     163         DO ji = mi0(ind1), mi1(ind2)    
     164            DO jj = 1, jpj                
     165               ztabramp(ji,jj) = 1._wp 
     166            END DO 
     167         END DO 
     168      ENDIF 
     169      IF( lk_east ) THEN                             ! --- East --- ! 
     170         ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1 
     171         ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     172         DO ji = mi0(ind1), mi1(ind2) 
     173            DO jj = 1, jpj 
     174               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     175            END DO 
     176         END DO 
     177         ! ghost cells: 
     178         ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     179         ind2 = jpiglo - 1 
     180         DO ji = mi0(ind1), mi1(ind2) 
     181            DO jj = 1, jpj 
     182               ztabramp(ji,jj) = 1._wp 
     183            END DO 
     184         END DO 
     185      ENDIF       
     186      IF( lk_south ) THEN                            ! --- South --- ! 
     187         ind1 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     188         ind2 = nn_hls + 1 + nbghostcells + jspongearea  
     189         DO jj = mj0(ind1), mj1(ind2)  
     190            DO ji = 1, jpi 
     191               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     192            END DO 
     193         END DO 
     194         ! ghost cells: 
     195         ind1 = 1 
     196         ind2 = nn_hls + 1 + nbghostcells                 ! halo + land + nbghostcells 
     197         DO jj = mj0(ind1), mj1(ind2)  
     198            DO ji = 1, jpi 
     199               ztabramp(ji,jj) = 1._wp 
     200            END DO 
     201         END DO 
     202      ENDIF 
     203      IF( lk_north ) THEN                            ! --- North --- ! 
     204         ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1 
     205         ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1    ! halo + land + nbghostcells - 1 
     206         DO jj = mj0(ind1), mj1(ind2) 
     207            DO ji = 1, jpi 
     208               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     209            END DO 
     210         END DO 
     211         ! ghost cells: 
     212         ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
     213         ind2 = jpjglo 
     214         DO jj = mj0(ind1), mj1(ind2) 
     215            DO ji = 1, jpi 
     216               ztabramp(ji,jj) = 1._wp 
     217            END DO 
     218         END DO 
     219      ENDIF       
     220      ! 
     221      ! Tracers 
     222      fspu(:,:) = 0._wp 
     223      fspv(:,:) = 0._wp 
     224      DO_2D( 0, 0, 0, 0 ) 
     225         fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     226         fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     227      END_2D 
     228 
     229      ! Dynamics 
     230      fspt(:,:) = 0._wp 
     231      fspf(:,:) = 0._wp 
     232      DO_2D( 0, 0, 0, 0 ) 
     233         fspt(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     234                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     235         fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
     236      END_2D 
     237       
     238      CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
     239      ! 
     240      ! Remove vertical interpolation where not needed: 
     241      ! (A null value in mbkx arrays does the job) 
     242      WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0 
     243      WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0 
     244      WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0 
     245      ! 
    182246#endif 
    183  
    184          ! Define ramp from boundaries towards domain interior at T-points 
    185          ! Store it in ztabramp 
    186  
    187          ispongearea    = nn_sponge_len * Agrif_irhox() 
    188          z1_ispongearea = 1._wp / REAL( ispongearea, wp ) 
    189          jspongearea    = nn_sponge_len * Agrif_irhoy() 
    190          z1_jspongearea = 1._wp / REAL( jspongearea, wp ) 
    191           
    192          ztabramp(:,:) = 0._wp 
    193  
    194          ! Trick to remove sponge in 2DV domains: 
    195          IF ( nbcellsx <= 3 ) ispongearea = -1 
    196          IF ( nbcellsy <= 3 ) jspongearea = -1 
    197  
    198          IF( lk_west ) THEN                             ! --- West --- ! 
    199             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    200             ind2 = nn_hls + 1 + nbghostcells + ispongearea  
    201             DO ji = mi0(ind1), mi1(ind2)    
    202                DO jj = 1, jpj                
    203                   ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea   * zmskwest(jj) 
    204                END DO 
    205             END DO 
    206             ! ghost cells: 
    207             ind1 = 1 
    208             ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    209             DO ji = mi0(ind1), mi1(ind2)    
    210                DO jj = 1, jpj                
    211                   ztabramp(ji,jj) = zmskwest(jj) 
    212                END DO 
    213             END DO 
    214          ENDIF 
    215          IF( lk_east ) THEN                             ! --- East --- ! 
    216             ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea 
    217             ind2 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    218             DO ji = mi0(ind1), mi1(ind2) 
    219                DO jj = 1, jpj 
    220                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) * zmskeast(jj) 
    221                END DO 
    222             END DO 
    223             ! ghost cells: 
    224             ind1 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    225             ind2 = jpiglo 
    226             DO ji = mi0(ind1), mi1(ind2) 
    227                DO jj = 1, jpj 
    228                   ztabramp(ji,jj) = zmskeast(jj) 
    229                END DO 
    230             END DO 
    231          ENDIF       
    232          IF( lk_south ) THEN                            ! --- South --- ! 
    233             ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    234             ind2 = nn_hls + 1 + nbghostcells + jspongearea  
    235             DO jj = mj0(ind1), mj1(ind2)  
    236                DO ji = 1, jpi 
    237                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) * zmsksouth(ji) 
    238                END DO 
    239             END DO 
    240             ! ghost cells: 
    241             ind1 = 1 
    242             ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells 
    243             DO jj = mj0(ind1), mj1(ind2)  
    244                DO ji = 1, jpi 
    245                   ztabramp(ji,jj) = zmsksouth(ji) 
    246                END DO 
    247             END DO 
    248          ENDIF 
    249          IF( lk_north ) THEN                            ! --- North --- ! 
    250             ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea 
    251             ind2 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    252             DO jj = mj0(ind1), mj1(ind2) 
    253                DO ji = 1, jpi 
    254                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) * zmsknorth(ji) 
    255                END DO 
    256             END DO 
    257             ! ghost cells: 
    258             ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1 
    259             ind2 = jpjglo 
    260             DO jj = mj0(ind1), mj1(ind2) 
    261                DO ji = 1, jpi 
    262                   ztabramp(ji,jj) = zmsknorth(ji) 
    263                END DO 
    264             END DO 
    265          ENDIF 
    266          ! 
     247      ! 
     248   END SUBROUTINE Agrif_Sponge 
     249 
     250 
     251   SUBROUTINE Agrif_Sponge_2d 
     252      !!---------------------------------------------------------------------- 
     253      !!                 *** ROUTINE  Agrif_Sponge_2d *** 
     254      !!---------------------------------------------------------------------- 
     255      INTEGER  ::   ji, jj, ind1, ind2, ishift, jshift 
     256      INTEGER  ::   ispongearea, jspongearea 
     257      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
     258      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
     259      !!---------------------------------------------------------------------- 
     260      ! 
     261      ! Sponge 1d example with: 
     262      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 
     263      !                         
     264      !coarse :     U     T     U     T     U     T     U 
     265      !|            |           |           |           | 
     266      !fine :     t u t u t u t u t u t u t u t u t u t u t 
     267      !sponge val:0 1   1   1   1  5/6 4/6 3/6 2/6 1/6  0  
     268      !           |   ghost     | <-- sponge area  -- > | 
     269      !           |   points    |                       | 
     270      !                         |--> dynamical interface 
     271 
     272#if defined SPONGE || defined SPONGE_TOP 
     273      ! Define ramp from boundaries towards domain interior at F-points 
     274      ! Store it in ztabramp 
     275 
     276      ispongearea    = nn_sponge_len * Agrif_irhox() 
     277      z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp ) 
     278      jspongearea    = nn_sponge_len * Agrif_irhoy() 
     279      z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp ) 
     280      ishift = nn_shift_bar * Agrif_irhox() 
     281      jshift = nn_shift_bar * Agrif_irhoy() 
     282         
     283      ztabramp(:,:) = 0._wp 
     284 
     285      IF( lk_west ) THEN                             ! --- West --- ! 
     286         ind1 = nn_hls + 1 + nbghostcells + ishift 
     287         ind2 = nn_hls + 1 + nbghostcells + ishift + ispongearea  
     288         DO ji = mi0(ind1), mi1(ind2)    
     289            DO jj = 1, jpj                
     290               ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea 
     291            END DO 
     292         END DO 
     293         ! ghost cells: 
     294         ind1 = 1 
     295         ind2 = nn_hls + 1 + nbghostcells + ishift               ! halo + land + nbghostcells 
     296         DO ji = mi0(ind1), mi1(ind2)    
     297            DO jj = 1, jpj                
     298               ztabramp(ji,jj) = 1._wp 
     299            END DO 
     300         END DO 
    267301      ENDIF 
    268  
     302      IF( lk_east ) THEN                             ! --- East --- ! 
     303         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - ispongearea - 1 
     304         ind2 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     305         DO ji = mi0(ind1), mi1(ind2) 
     306            DO jj = 1, jpj 
     307               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )  
     308            END DO 
     309         END DO 
     310         ! ghost cells: 
     311         ind1 = jpiglo - ( nn_hls + nbghostcells + ishift) - 1    ! halo + land + nbghostcells - 1 
     312         ind2 = jpiglo - 1 
     313         DO ji = mi0(ind1), mi1(ind2) 
     314            DO jj = 1, jpj 
     315               ztabramp(ji,jj) = 1._wp 
     316            END DO 
     317         END DO 
     318      ENDIF       
     319      IF( lk_south ) THEN                            ! --- South --- ! 
     320         ind1 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     321         ind2 = nn_hls + 1 + nbghostcells + jshift + jspongearea  
     322         DO jj = mj0(ind1), mj1(ind2)  
     323            DO ji = 1, jpi 
     324               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) 
     325            END DO 
     326         END DO 
     327         ! ghost cells: 
     328         ind1 = 1 
     329         ind2 = nn_hls + 1 + nbghostcells + jshift                ! halo + land + nbghostcells 
     330         DO jj = mj0(ind1), mj1(ind2)  
     331            DO ji = 1, jpi 
     332               ztabramp(ji,jj) = 1._wp 
     333            END DO 
     334         END DO 
     335      ENDIF 
     336      IF( lk_north ) THEN                            ! --- North --- ! 
     337         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift) - jspongearea - 1 
     338         ind2 = jpjglo - ( nn_hls + nbghostcells + jshift) - 1    ! halo + land + nbghostcells - 1 
     339         DO jj = mj0(ind1), mj1(ind2) 
     340            DO ji = 1, jpi 
     341               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )  
     342            END DO 
     343         END DO 
     344         ! ghost cells: 
     345         ind1 = jpjglo - ( nn_hls + nbghostcells + jshift)      ! halo + land + nbghostcells - 1 
     346         ind2 = jpjglo 
     347         DO jj = mj0(ind1), mj1(ind2) 
     348            DO ji = 1, jpi 
     349               ztabramp(ji,jj) = 1._wp 
     350            END DO 
     351         END DO 
     352      ENDIF 
     353      ! 
    269354      ! Tracers 
    270       IF( .NOT. spongedoneT ) THEN 
    271          fspu(:,:) = 0._wp 
    272          fspv(:,:) = 0._wp 
    273          DO_2D( 0, 0, 0, 0 ) 
    274             fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) * ssumask(ji,jj) 
    275             fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) * ssvmask(ji,jj) 
     355      fspu_2d(:,:) = 0._wp 
     356      fspv_2d(:,:) = 0._wp 
     357      DO_2D( 0, 0, 0, 0 ) 
     358         fspu_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj) 
     359         fspv_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj) 
     360      END_2D 
     361 
     362      ! Dynamics 
     363      fspt_2d(:,:) = 0._wp 
     364      fspf_2d(:,:) = 0._wp 
     365      DO_2D( 0, 0, 0, 0 ) 
     366         fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji-1,jj  ) & 
     367                               &  +ztabramp(ji  ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj) 
     368         fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
    276369         END_2D 
    277       ENDIF 
    278  
    279       ! Dynamics 
    280       IF( .NOT. spongedoneU ) THEN 
    281          fspt(:,:) = 0._wp 
    282          fspf(:,:) = 0._wp 
    283          DO_2D( 0, 0, 0, 0 ) 
    284             fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 
    285             fspf(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1)   & 
    286                                   &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) ) & 
    287                                   &  * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
    288          END_2D 
    289       ENDIF 
    290        
    291       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    292          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    293          spongedoneT = .TRUE. 
    294          spongedoneU = .TRUE. 
    295       ENDIF 
    296       IF( .NOT. spongedoneT ) THEN 
    297          CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp ) 
    298          spongedoneT = .TRUE. 
    299       ENDIF 
    300       IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 
    301          CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp ) 
    302          spongedoneU = .TRUE. 
    303       ENDIF 
    304  
    305 #if defined key_vertical 
    306       ! Remove vertical interpolation where not needed: 
    307       DO_2D( 0, 0, 0, 0 ) 
    308          IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 
    309          &   (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 
    310 ! 
    311          IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
    312          &   (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 
    313 ! 
    314          IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
    315          &   (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 
    316 ! 
    317          IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 
    318          IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 
    319          IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 
    320       END_2D 
    321       ! 
    322       ztabramp (:,:) = REAL( mbkt_parent(:,:), wp ) 
    323       ztabrampu(:,:) = REAL( mbku_parent(:,:), wp ) 
    324       ztabrampv(:,:) = REAL( mbkv_parent(:,:), wp ) 
    325       CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1._wp, ztabrampu, 'U', 1._wp, ztabrampv, 'V', 1._wp ) 
    326       mbkt_parent(:,:) = NINT( ztabramp (:,:) ) 
    327       mbku_parent(:,:) = NINT( ztabrampu(:,:) ) 
    328       mbkv_parent(:,:) = NINT( ztabrampv(:,:) ) 
     370      CALL lbc_lnk( 'agrif_Sponge_2d', fspu_2d, 'U', 1._wp, fspv_2d, 'V', 1._wp, fspt_2d, 'T', 1._wp, fspf_2d, 'F', 1._wp ) 
     371      ! 
    329372#endif 
    330373      ! 
    331 #endif 
    332       ! 
    333    END SUBROUTINE Agrif_Sponge 
    334  
    335     
    336    SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     374   END SUBROUTINE Agrif_Sponge_2d 
     375 
     376 
     377   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)  
    337378      !!---------------------------------------------------------------------- 
    338379      !!                 *** ROUTINE interptsn_sponge *** 
    … …  
    345386      INTEGER  ::   iku, ikv 
    346387      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 
    347       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
     388      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv 
    348389      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
    349390      ! vertical interpolation: 
    350391      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
    351       REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
    352       REAL(wp), DIMENSION(k1:k2) :: h_in 
    353       REAL(wp), DIMENSION(1:jpk) :: h_out 
     392      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i 
     393      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i 
     394      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    354395      INTEGER :: N_in, N_out 
    355396      !!---------------------------------------------------------------------- 
    … …  
    366407         END DO 
    367408 
    368 # if defined key_vertical 
    369         ! Interpolate thicknesses 
    370         ! Warning: these are masked, hence extrapolated prior interpolation. 
    371         DO jk=k1,k2 
    372            DO jj=j1,j2 
    373               DO ji=i1,i2 
    374                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 
    375               END DO 
    376            END DO 
    377         END DO 
    378  
    379         ! Extrapolate thicknesses in partial bottom cells: 
    380         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    381         IF (ln_zps) THEN 
    382            DO jj=j1,j2 
    383               DO ji=i1,i2 
    384                   jk = mbkt(ji,jj) 
    385                   tabres(ji,jj,jk,jpts+1) = 0._wp 
    386               END DO 
    387            END DO            
    388         END IF 
    389       
    390         ! Save ssh at last level: 
    391         IF (.NOT.ln_linssh) THEN 
    392            tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
    393         ELSE 
    394            tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
    395         END IF       
    396 # endif 
    397  
    398       ELSE    
    399          ! 
    400 # if defined key_vertical 
    401  
    402          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
    403  
    404          DO jj=j1,j2 
    405             DO ji=i1,i2 
    406                tabres_child(ji,jj,:,:) = 0._wp  
    407                N_in = mbkt_parent(ji,jj) 
    408                zhtot = 0._wp 
    409                DO jk=1,N_in !k2 = jpk of parent grid 
    410                   IF (jk==N_in) THEN 
    411                      h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 
     409         IF ( l_vremap.OR.ln_zps ) THEN 
     410 
     411            ! Fill cell depths (i.e. gdept) to be interpolated 
     412            ! Warning: these are masked, hence extrapolated prior interpolation. 
     413            DO jj=j1,j2 
     414               DO ji=i1,i2 
     415                  tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a) 
     416                  DO jk=k1+1,k2 
     417                     tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * & 
     418                        & ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) ) 
     419                  END DO 
     420               END DO 
     421            END DO 
     422 
     423            ! Save ssh at last level: 
     424            IF ( .NOT.ln_linssh ) THEN 
     425               tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
     426            END IF   
     427     
     428         END IF 
     429 
     430      ELSE  
     431         ! 
     432         IF ( l_vremap ) THEN 
     433 
     434            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
     435 
     436            DO jj=j1,j2 
     437               DO ji=i1,i2 
     438 
     439                  tabres_child(ji,jj,:,:) = 0._wp  
     440                  ! Build vertical grids: 
     441                  N_in = mbkt_parent(ji,jj) 
     442                  ! Input grid (account for partial cells if any): 
     443                  IF ( N_in > 0 ) THEN  
     444                     DO jk=1,N_in 
     445                        z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2) 
     446                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     447                     END DO 
     448                   
     449                     ! Intermediate grid: 
     450                     DO jk = 1, N_in 
     451                        h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     452                          &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     453                     END DO 
     454                     z_in_i(1) = 0.5_wp * h_in_i(1) 
     455                     DO jk=2,N_in 
     456                        z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     457                     END DO 
     458                     z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2) 
     459                  END IF 
     460                  ! Output (Child) grid: 
     461                  N_out = mbkt(ji,jj) 
     462                  DO jk=1,N_out 
     463                     h_out(jk) = e3t(ji,jj,jk,Kbb_a) 
     464                  END DO 
     465                  z_out(1) = 0.5_wp * h_out(1) 
     466                  DO jk=2,N_out 
     467                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     468                  END DO 
     469                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a) 
     470 
     471                  ! Account for small differences in the free-surface 
     472                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN 
     473                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) ) 
    412474                  ELSE 
    413                      h_in(jk) = tabres(ji,jj,jk,n2) 
     475                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) ) 
     476                  END IF 
     477                  IF (N_in*N_out > 0) THEN 
     478                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),z_in_i(1:N_in),N_in,N_in,jpts) 
     479                     CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
     480!                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),z_out(1:N_in),N_in,N_out,jpts)   
    414481                  ENDIF 
    415                   zhtot = zhtot + h_in(jk) 
    416                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    417                END DO 
    418                N_out = 0 
    419                DO jk=1,jpk ! jpk of child grid 
    420                   IF (tmask(ji,jj,jk) == 0) EXIT  
    421                   N_out = N_out + 1 
    422                   h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    423                END DO 
    424  
    425                ! Account for small differences in free-surface 
    426                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    427                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    428                ELSE 
    429                   h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    430                ENDIF 
    431                IF (N_in*N_out > 0) THEN 
    432                   CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    433                ENDIF 
    434             END DO 
    435          END DO 
    436 # endif 
    437  
    438          DO jj=j1,j2 
    439             DO ji=i1,i2 
    440                DO jk=1,jpkm1 
    441 # if defined key_vertical 
    442                   tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 
    443 # else 
    444                   tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 
    445 # endif 
    446                END DO 
    447             END DO 
    448          END DO 
     482               END DO 
     483            END DO 
     484 
     485            DO jj=j1,j2 
     486               DO ji=i1,i2 
     487                  DO jk=1,jpkm1 
     488                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 
     489                  END DO 
     490               END DO 
     491            END DO 
     492 
     493         ELSE 
     494 
     495            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
     496 
     497               DO jj=j1,j2 
     498                  DO ji=i1,i2 
     499                     ! 
     500                     N_in  = mbkt(ji,jj)  
     501                     N_out = mbkt(ji,jj)  
     502                     z_in(1) = tabres(ji,jj,1,n2) 
     503                     tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts) 
     504                     DO jk=2, N_in 
     505                        z_in(jk) = tabres(ji,jj,jk,n2) 
     506                        tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts) 
     507                     END DO  
     508                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2) 
     509 
     510                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a) 
     511                     DO jk=2, N_out 
     512                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a))  
     513                     END DO  
     514                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a) 
     515 
     516                     CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), & 
     517                                         &   z_out(1:N_out), N_in, N_out, jpts) 
     518                  END DO 
     519               END DO 
     520            ENDIF 
     521 
     522            DO jj=j1,j2 
     523               DO ji=i1,i2 
     524                  DO jk=1,jpkm1 
     525                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 
     526                  END DO 
     527               END DO 
     528            END DO 
     529 
     530         END IF 
    449531 
    450532         DO jn = 1, jpts             
    451533            DO jk = 1, jpkm1 
    452                ztu(i1:i2,j1:j2,jk) = 0._wp 
     534               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp 
    453535               DO jj = j1,j2 
    454536                  DO ji = i1,i2-1 
    455                      zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    456                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
    457                   END DO 
    458                END DO 
    459                ztv(i1:i2,j1:j2,jk) = 0._wp 
     537                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
     538                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
     539                  END DO 
     540               END DO 
     541               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp 
    460542               DO ji = i1,i2 
    461543                  DO jj = j1,j2-1 
    462                      zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    463                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     544                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
     545                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    464546                  END DO 
    465547               END DO 
    … …  
    478560            END DO 
    479561            ! 
     562! JC: there is something wrong with the Laplacian in corners 
    480563            DO jk = 1, jpkm1 
    481                DO jj = j1+1,j2-1 
    482                   DO ji = i1+1,i2-1 
     564               DO jj = j1,j2 
     565                  DO ji = i1,i2 
    483566                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN  
    484567                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 
    485568                        ! horizontal diffusive trends 
    486                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) & 
    487                              &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)  
     569                        ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &  
     570                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 
     571                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
    488572                        ! add it to the general tracer trends 
    489573                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa 
    … …  
    491575                  END DO 
    492576               END DO 
     577 
    493578            END DO 
    494579            ! 
    495580         END DO 
    496581         ! 
    497          tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     582         tabspongedone_tsn(i1:i2,j1:j2) = .TRUE. 
    498583         ! 
    499584      ENDIF 
    … …  
    528613               DO ji=i1,i2 
    529614                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 
    530 # if defined key_vertical 
    531                   tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 
    532 # endif 
    533                END DO 
    534             END DO 
    535          END DO 
    536  
    537 # if defined key_vertical 
    538          ! Extrapolate thicknesses in partial bottom cells: 
    539          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    540          IF (ln_zps) THEN 
     615               END DO 
     616            END DO 
     617         END DO 
     618 
     619         IF ( l_vremap ) THEN 
     620 
     621            DO jk=k1,k2 
     622               DO jj=j1,j2 
     623                  DO ji=i1,i2 
     624                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 
     625                  END DO 
     626               END DO 
     627            END DO 
     628 
     629            ! Extrapolate thicknesses in partial bottom cells: 
     630            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     631            IF (ln_zps) THEN 
     632               DO jj=j1,j2 
     633                  DO ji=i1,i2 
     634                     jk = mbku(ji,jj) 
     635                     tabres(ji,jj,jk,m2) = 0._wp 
     636                  END DO 
     637               END DO            
     638            END IF 
     639            ! Save ssh at last level: 
     640            tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     641            IF (.NOT.ln_linssh) THEN 
     642               ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
     643               DO jk=1,jpk 
     644                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 
     645               END DO 
     646               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 
     647            END IF  
     648         END IF 
     649 
     650      ELSE 
     651 
     652         IF ( l_vremap ) THEN 
     653 
     654            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     655 
    541656            DO jj=j1,j2 
    542657               DO ji=i1,i2 
    543                   jk = mbku(ji,jj) 
    544                   tabres(ji,jj,jk,m2) = 0._wp 
    545                END DO 
    546             END DO            
    547          END IF 
    548         ! Save ssh at last level: 
    549         tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    550         IF (.NOT.ln_linssh) THEN 
    551            ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
    552            DO jk=1,jpk 
    553               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 
    554            END DO 
    555            tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 
    556         END IF  
    557 # endif 
    558  
    559       ELSE 
    560  
    561 # if defined key_vertical 
    562          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    563  
    564          DO jj=j1,j2 
    565             DO ji=i1,i2 
    566                tabres_child(ji,jj,:) = 0._wp 
    567                N_in = mbku_parent(ji,jj) 
    568                zhtot = 0._wp 
    569                DO jk=1,N_in 
    570                   IF (jk==N_in) THEN 
    571                      h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
    572                   ELSE 
    573                      h_in(jk) = tabres(ji,jj,jk,m2) 
    574                   ENDIF 
    575                   zhtot = zhtot + h_in(jk) 
    576                   tabin(jk) = tabres(ji,jj,jk,m1) 
    577                END DO 
    578                !          
    579                N_out = 0 
    580                DO jk=1,jpk 
    581                   IF (umask(ji,jj,jk) == 0) EXIT 
    582                   N_out = N_out + 1 
    583                   h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 
    584                END DO 
    585  
    586                ! Account for small differences in free-surface 
    587                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    588                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    589                ELSE 
    590                   h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    591                ENDIF 
     658                  tabres_child(ji,jj,:) = 0._wp 
     659                  N_in = mbku_parent(ji,jj) 
     660                  N_out = mbku(ji,jj) 
     661                  IF (N_in * N_out > 0) THEN 
     662                     zhtot = 0._wp 
     663                     DO jk=1,N_in 
     664                        !IF (jk==N_in) THEN 
     665                        !   h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     666                        !ELSE 
     667                        !   h_in(jk) = tabres(ji,jj,jk,m2) 
     668                        !ENDIF 
     669                        h_in(jk) = e3u0_parent(ji,jj,jk) 
     670                        zhtot = zhtot + h_in(jk) 
     671                        tabin(jk) = tabres(ji,jj,jk,m1) 
     672                     END DO 
     673                     !          
     674                     DO jk=1,N_out 
     675                        h_out(jk) = e3u(ji,jj,jk,Kbb_a) 
     676                     END DO 
     677 
     678                     ! Account for small differences in free-surface 
     679                     IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     680                        h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     681                     ELSE 
     682                        h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     683                     ENDIF 
    592684                   
    593                IF (N_in * N_out > 0) THEN 
    594                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    595                ENDIF  
    596             END DO 
    597          END DO 
    598  
    599          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
    600 #else 
    601          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
    602 #endif 
     685                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
     686                  ENDIF  
     687               END DO 
     688            END DO 
     689 
     690            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
     691         ELSE 
     692 
     693            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
     694   
     695         ENDIF 
    603696         ! 
    604697         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    705798               DO ji=i1,i2 
    706799                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) 
    707 # if defined key_vertical 
    708                   tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 
    709 # endif 
    710                END DO 
    711             END DO 
    712          END DO 
    713  
    714 # if defined key_vertical 
    715          ! Extrapolate thicknesses in partial bottom cells: 
    716          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    717          IF (ln_zps) THEN 
     800               END DO 
     801            END DO 
     802         END DO 
     803 
     804         IF ( l_vremap ) THEN 
     805 
     806            DO jk=k1,k2 
     807               DO jj=j1,j2 
     808                  DO ji=i1,i2 
     809                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 
     810                  END DO 
     811               END DO 
     812            END DO 
     813            ! Extrapolate thicknesses in partial bottom cells: 
     814            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     815            IF (ln_zps) THEN 
     816               DO jj=j1,j2 
     817                  DO ji=i1,i2 
     818                     jk = mbkv(ji,jj) 
     819                     tabres(ji,jj,jk,m2) = 0._wp 
     820                  END DO 
     821               END DO            
     822            END IF 
     823            ! Save ssh at last level: 
     824            tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     825            IF (.NOT.ln_linssh) THEN 
     826               ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
     827               DO jk=1,jpk 
     828                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk) 
     829               END DO 
     830               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 
     831            END IF 
     832 
     833         END IF  
     834 
     835      ELSE 
     836 
     837         IF ( l_vremap ) THEN 
     838            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    718839            DO jj=j1,j2 
    719840               DO ji=i1,i2 
    720                   jk = mbkv(ji,jj) 
    721                   tabres(ji,jj,jk,m2) = 0._wp 
    722                END DO 
    723             END DO            
    724          END IF 
    725         ! Save ssh at last level: 
    726         tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    727         IF (.NOT.ln_linssh) THEN 
    728            ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
    729            DO jk=1,jpk 
    730               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk) 
    731            END DO 
    732            tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 
    733         END IF  
    734 # endif 
    735  
    736       ELSE 
    737  
    738 # if defined key_vertical 
    739          IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    740          DO jj=j1,j2 
    741             DO ji=i1,i2 
    742                tabres_child(ji,jj,:) = 0._wp 
    743                N_in = mbkv_parent(ji,jj) 
    744                zhtot = 0._wp 
    745                DO jk=1,N_in 
    746                   IF (jk==N_in) THEN 
    747                      h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
    748                   ELSE 
    749                      h_in(jk) = tabres(ji,jj,jk,m2) 
     841                  tabres_child(ji,jj,:) = 0._wp 
     842                  N_in = mbkv_parent(ji,jj) 
     843                  N_out = mbkv(ji,jj) 
     844                  IF (N_in * N_out > 0) THEN 
     845                     zhtot = 0._wp 
     846                     DO jk=1,N_in 
     847                        !IF (jk==N_in) THEN 
     848                        !   h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     849                        !ELSE 
     850                        !   h_in(jk) = tabres(ji,jj,jk,m2) 
     851                        !ENDIF 
     852                        h_in(jk) = e3v0_parent(ji,jj,jk) 
     853                        zhtot = zhtot + h_in(jk) 
     854                        tabin(jk) = tabres(ji,jj,jk,m1) 
     855                     END DO 
     856                     !           
     857                     DO jk=1,N_out 
     858                        h_out(jk) = e3v(ji,jj,jk,Kbb_a) 
     859                     END DO 
     860 
     861                     ! Account for small differences in free-surface 
     862                     IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     863                        h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     864                     ELSE 
     865                        h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     866                     ENDIF 
     867          
     868                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
     869 
    750870                  ENDIF 
    751                   zhtot = zhtot + h_in(jk) 
    752                   tabin(jk) = tabres(ji,jj,jk,m1) 
    753                END DO 
    754                !           
    755                N_out = 0 
    756                DO jk=1,jpk 
    757                   IF (vmask(ji,jj,jk) == 0) EXIT 
    758                   N_out = N_out + 1 
    759                   h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 
    760                END DO 
    761  
    762                ! Account for small differences in free-surface 
    763                IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
    764                   h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
    765                ELSE 
    766                   h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
    767                ENDIF 
    768           
    769                IF (N_in * N_out > 0) THEN 
    770                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    771                ENDIF 
    772             END DO 
    773          END DO 
    774  
    775          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
    776 # else 
    777          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
    778 # endif 
     871               END DO 
     872            END DO 
     873 
     874            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
     875         ELSE 
     876 
     877            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
     878 
     879         ENDIF 
    779880         ! 
    780881         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    839940      ! 
    840941   END SUBROUTINE interpvn_sponge 
     942 
     943   SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before) 
     944      !!--------------------------------------------- 
     945      !!   *** ROUTINE interpunb_sponge *** 
     946      !!---------------------------------------------     
     947      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     948      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     949      LOGICAL, INTENT(in) :: before 
     950 
     951      INTEGER  :: ji, jj, ind1, jmax 
     952      ! sponge parameters  
     953      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     954      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff 
     955      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     956      !!---------------------------------------------     
     957      ! 
     958      IF( before ) THEN 
     959         DO jj=j1,j2 
     960            DO ji=i1,i2 
     961               tabres(ji,jj) = uu_b(ji,jj,Kmm_a) 
     962            END DO 
     963         END DO 
     964 
     965      ELSE 
     966 
     967         ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1) 
     968         ! 
     969         !                                             ! -------- 
     970         ! Horizontal divergence                       !   div 
     971         !                                             ! -------- 
     972         DO jj = j1,j2 
     973            DO ji = i1+1,i2   ! vector opt. 
     974               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     975               hdivdiff(ji,jj) = (  e2u(ji  ,jj)*hu(ji  ,jj,Kbb_a) * ubdiff(ji  ,jj) & 
     976                                  &-e2u(ji-1,jj)*hu(ji-1,jj,Kbb_a) * ubdiff(ji-1,jj) ) * zbtr 
     977            END DO 
     978         END DO 
     979 
     980         DO jj = j1,j2-1 
     981            DO ji = i1,i2   ! vector opt. 
     982               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj) 
     983               rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1)   & 
     984                              &   +e1u(ji,jj  ) * ubdiff(ji,jj  ) ) * fmask(ji,jj,1) * zbtr  
     985            END DO 
     986         END DO 
     987         ! 
     988         DO jj = j1+1, j2-1 
     989            DO ji = i1+1, i2-1   ! vector opt. 
     990               IF (.NOT. tabspongedone_u(ji,jj)) THEN 
     991                  ze2u = rotdiff (ji,jj) 
     992                  ze1v = hdivdiff(ji,jj) 
     993                  ! horizontal diffusive trends 
     994                  zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     995                      & + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj)                       &  
     996                      & - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj) 
     997 
     998                  ! add it to the general momentum trends 
     999                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua                                  
     1000               ENDIF 
     1001            END DO 
     1002         END DO 
     1003 
     1004         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1005 
     1006         jmax = j2-1 
     1007         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North 
     1008         DO jj = mj0(ind1), mj1(ind1)                  
     1009            jmax = MIN(jmax,jj) 
     1010         END DO 
     1011 
     1012         DO jj = j1+1, jmax 
     1013            DO ji = i1+1, i2   ! vector opt. 
     1014               IF (.NOT. tabspongedone_v(ji,jj)) THEN 
     1015                     ze2u = rotdiff (ji,jj) 
     1016                     ze1v = hdivdiff(ji,jj) 
     1017                     zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) & 
     1018                           + ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj) 
     1019                     vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1020               ENDIF 
     1021            END DO 
     1022         END DO 
     1023         ! 
     1024         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE. 
     1025         ! 
     1026      ENDIF 
     1027      ! 
     1028   END SUBROUTINE interpunb_sponge 
     1029 
     1030    
     1031   SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before) 
     1032      !!--------------------------------------------- 
     1033      !!   *** ROUTINE interpvnb_sponge *** 
     1034      !!---------------------------------------------  
     1035      INTEGER, INTENT(in) :: i1,i2,j1,j2 
     1036      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1037      LOGICAL, INTENT(in) :: before 
     1038      ! 
     1039      INTEGER  ::   ji, jj, ind1, imax 
     1040      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr 
     1041      REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff 
     1042      REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff 
     1043      !!---------------------------------------------  
     1044       
     1045      IF( before ) THEN  
     1046         DO jj=j1,j2 
     1047            DO ji=i1,i2 
     1048               tabres(ji,jj) = vv_b(ji,jj,Kmm_a) 
     1049            END DO 
     1050         END DO 
     1051      ELSE 
     1052         vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1) 
     1053         !                                             ! -------- 
     1054         ! Horizontal divergence                       !   div 
     1055         !                                             ! -------- 
     1056         DO jj = j1+1,j2 
     1057            DO ji = i1,i2   ! vector opt. 
     1058               zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj) 
     1059               hdivdiff(ji,jj) = ( e1v(ji,jj  ) * hv(ji,jj  ,Kbb_a) * vbdiff(ji,jj  )  & 
     1060                               &  -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1)  ) * zbtr 
     1061            END DO 
     1062         END DO 
     1063         DO jj = j1,j2 
     1064            DO ji = i1,i2-1   ! vector opt. 
     1065               zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)  
     1066               rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) &  
     1067                              &  -e2v(ji  ,jj) * vbdiff(ji  ,jj)  ) * fmask(ji,jj,1) * zbtr 
     1068            END DO 
     1069         END DO 
     1070         !                                                ! =============== 
     1071         !                                                 
     1072 
     1073         imax = i2 - 1 
     1074         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East 
     1075         DO ji = mi0(ind1), mi1(ind1)                 
     1076            imax = MIN(imax,ji) 
     1077         END DO 
     1078          
     1079         DO jj = j1+1, j2 
     1080            DO ji = i1+1, imax   ! vector opt. 
     1081               IF( .NOT. tabspongedone_u(ji,jj) ) THEN                                                      
     1082                  zua = - ( rotdiff (ji  ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a)  & 
     1083                      & + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj  )) * r1_e1u(ji,jj) 
     1084                  uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua 
     1085               ENDIF 
     1086            END DO 
     1087         END DO 
     1088         ! 
     1089         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE. 
     1090         ! 
     1091         DO jj = j1+1, j2-1 
     1092            DO ji = i1+1, i2-1   ! vector opt. 
     1093               IF( .NOT. tabspongedone_v(ji,jj) ) THEN 
     1094                  zva  =  ( rotdiff (ji,jj  ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) & 
     1095                     &  + ( hdivdiff(ji,jj+1) - hdivdiff(ji  ,jj) ) * r1_e2v(ji,jj)                     & 
     1096                     &  - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj) 
     1097                  vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva 
     1098               ENDIF 
     1099            END DO 
     1100         END DO 
     1101         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     1102      ENDIF 
     1103      ! 
     1104   END SUBROUTINE interpvnb_sponge 
     1105 
    8411106 
    8421107#else 
Note: See TracChangeset for help on using the changeset viewer.