New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12377 for NEMO/trunk/src/NST/agrif_oce_sponge.F90 – NEMO

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

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

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

svn merge --ignore-ancestry \

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

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

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

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

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

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

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

    r10425 r12377  
    2222   USE agrif_oce 
    2323   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     24   USE iom 
     25   USE vremap 
    2426 
    2527   IMPLICIT NONE 
    … …  
    2931   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge 
    3032 
     33   !! * Substitutions 
     34#  include "do_loop_substitute.h90" 
    3135   !!---------------------------------------------------------------------- 
    3236   !! NEMO/NST 4.0 , NEMO Consortium (2018) 
    … …  
    5862#endif 
    5963      ! 
     64      CALL iom_put( 'agrif_spu', fspu(:,:)) 
     65      CALL iom_put( 'agrif_spv', fspv(:,:)) 
     66      ! 
    6067   END SUBROUTINE Agrif_Sponge_Tra 
    6168 
    … …  
    8592#endif 
    8693      ! 
     94      CALL iom_put( 'agrif_spt', fspt(:,:)) 
     95      CALL iom_put( 'agrif_spf', fspf(:,:)) 
     96      ! 
    8797   END SUBROUTINE Agrif_Sponge_dyn 
    8898 
    … …  
    93103      !!---------------------------------------------------------------------- 
    94104      INTEGER  ::   ji, jj, ind1, ind2 
    95       INTEGER  ::   ispongearea 
    96       REAL(wp) ::   z1_spongearea 
     105      INTEGER  ::   ispongearea, jspongearea 
     106      REAL(wp) ::   z1_ispongearea, z1_jspongearea 
    97107      REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 
    98       !!---------------------------------------------------------------------- 
    99       ! 
     108      REAL(wp), DIMENSION(jpjmax)  :: zmskwest,  zmskeast 
     109      REAL(wp), DIMENSION(jpimax)  :: zmsknorth, zmsksouth 
     110      !!---------------------------------------------------------------------- 
     111      ! 
     112      ! Sponge 1d example with: 
     113      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 
     114      !                         
     115      !coarse :     U     T     U     T     U     T     U 
     116      !|            |           |           |           | 
     117      !fine :     t u t u t u t u t u t u t u t u t u t u t 
     118      !sponge val:0   0   0   1  5/6 4/6 3/6 2/6 1/6  0   0 
     119      !           |   ghost     | <-- sponge area  -- > | 
     120      !           |   points    |                       | 
     121      !                         |--> dynamical interface 
     122 
    100123#if defined SPONGE || defined SPONGE_TOP 
    101124      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
     125         ! 
     126         ! Retrieve masks at open boundaries: 
     127 
     128         ! --- West --- ! 
     129         ztabramp(:,:) = 0._wp 
     130         ind1 = 1+nbghostcells 
     131         DO ji = mi0(ind1), mi1(ind1)                 
     132            ztabramp(ji,:) = ssumask(ji,:) 
     133         END DO 
     134         ! 
     135         zmskwest(:) = 0._wp 
     136         zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 
     137 
     138         ! --- East --- ! 
     139         ztabramp(:,:) = 0._wp 
     140         ind1 = jpiglo - nbghostcells - 1 
     141         DO ji = mi0(ind1), mi1(ind1)                  
     142            ztabramp(ji,:) = ssumask(ji,:) 
     143         END DO 
     144         ! 
     145         zmskeast(:) = 0._wp 
     146         zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 
     147 
     148         ! --- South --- ! 
     149         ztabramp(:,:) = 0._wp 
     150         ind1 = 1+nbghostcells 
     151         DO jj = mj0(ind1), mj1(ind1)                  
     152            ztabramp(:,jj) = ssvmask(:,jj) 
     153         END DO 
     154         ! 
     155         zmsksouth(:) = 0._wp 
     156         zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 
     157 
     158         ! --- North --- ! 
     159         ztabramp(:,:) = 0._wp 
     160         ind1 = jpjglo - nbghostcells - 1 
     161         DO jj = mj0(ind1), mj1(ind1)                  
     162            ztabramp(:,jj) = ssvmask(:,jj) 
     163         END DO 
     164         ! 
     165         zmsknorth(:) = 0._wp 
     166         zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 
     167         ! JC: SPONGE MASKING TO BE SORTED OUT: 
     168         zmskwest(:)  = 1._wp 
     169         zmskeast(:)  = 1._wp 
     170         zmsknorth(:) = 1._wp 
     171         zmsksouth(:) = 1._wp 
     172#if defined key_mpp_mpi 
     173!         CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 
     174!         CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 
     175!         CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 
     176!         CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 
     177#endif 
     178 
    102179         ! Define ramp from boundaries towards domain interior at T-points 
    103180         ! Store it in ztabramp 
    104181 
    105          ispongearea  = 1 + nn_sponge_len * Agrif_irhox() 
    106          z1_spongearea = 1._wp / REAL( ispongearea ) 
     182         ispongearea  = nn_sponge_len * Agrif_irhox() 
     183         z1_ispongearea = 1._wp / REAL( ispongearea ) 
     184         jspongearea  = nn_sponge_len * Agrif_irhoy() 
     185         z1_jspongearea = 1._wp / REAL( jspongearea ) 
    107186          
    108187         ztabramp(:,:) = 0._wp 
    109188 
     189         ! Trick to remove sponge in 2DV domains: 
     190         IF ( nbcellsx <= 3 ) ispongearea = -1 
     191         IF ( nbcellsy <= 3 ) jspongearea = -1 
     192 
    110193         ! --- West --- ! 
    111          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    112             ind1 = 1+nbghostcells 
    113             ind2 = 1+nbghostcells + ispongearea  
     194         ind1 = 1+nbghostcells 
     195         ind2 = 1+nbghostcells + ispongearea  
     196         DO ji = mi0(ind1), mi1(ind2)    
     197            DO jj = 1, jpj                
     198               ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 
     199            END DO 
     200         END DO 
     201 
     202         ! ghost cells: 
     203         ind1 = 1 
     204         ind2 = nbghostcells + 1 
     205         DO ji = mi0(ind1), mi1(ind2)    
     206            DO jj = 1, jpj                
     207               ztabramp(ji,jj) = zmskwest(jj) 
     208            END DO 
     209         END DO 
     210 
     211         ! --- East --- ! 
     212         ind1 = jpiglo - nbghostcells - ispongearea 
     213         ind2 = jpiglo - nbghostcells 
     214         DO ji = mi0(ind1), mi1(ind2) 
    114215            DO jj = 1, jpj 
    115                DO ji = ind1, ind2                 
    116                   ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1) 
    117                END DO 
     216               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 
    118217            ENDDO 
    119          ENDIF 
    120  
    121          ! --- East --- ! 
    122          IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    123             ind1 = nlci - nbghostcells - ispongearea 
    124             ind2 = nlci - nbghostcells 
     218         END DO 
     219 
     220         ! ghost cells: 
     221         ind1 = jpiglo - nbghostcells 
     222         ind2 = jpiglo 
     223         DO ji = mi0(ind1), mi1(ind2) 
    125224            DO jj = 1, jpj 
    126                DO ji = ind1, ind2 
    127                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind1 ) * z1_spongearea * umask(ind2-1,jj,1) ) 
    128                ENDDO 
     225               ztabramp(ji,jj) = zmskeast(jj) 
    129226            ENDDO 
    130          ENDIF 
     227         END DO 
    131228 
    132229         ! --- South --- ! 
    133          IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    134             ind1 = 1+nbghostcells 
    135             ind2 = 1+nbghostcells + ispongearea 
    136             DO jj = ind1, ind2  
    137                DO ji = 1, jpi 
    138                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) ) 
    139                END DO 
    140             ENDDO 
    141          ENDIF 
     230         ind1 = 1+nbghostcells 
     231         ind2 = 1+nbghostcells + jspongearea 
     232         DO jj = mj0(ind1), mj1(ind2)  
     233            DO ji = 1, jpi 
     234               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 
     235            END DO 
     236         END DO 
     237 
     238         ! ghost cells: 
     239         ind1 = 1 
     240         ind2 = nbghostcells + 1 
     241         DO jj = mj0(ind1), mj1(ind2)  
     242            DO ji = 1, jpi 
     243               ztabramp(ji,jj) = zmsksouth(ji) 
     244            END DO 
     245         END DO 
    142246 
    143247         ! --- North --- ! 
    144          IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    145             ind1 = nlcj - nbghostcells - ispongearea 
    146             ind2 = nlcj - nbghostcells 
    147             DO jj = ind1, ind2 
    148                DO ji = 1, jpi 
    149                   ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind1 ) * z1_spongearea * vmask(ji,ind2-1,1) ) 
    150                END DO 
    151             ENDDO 
    152          ENDIF 
     248         ind1 = jpjglo - nbghostcells - jspongearea 
     249         ind2 = jpjglo - nbghostcells 
     250         DO jj = mj0(ind1), mj1(ind2) 
     251            DO ji = 1, jpi 
     252               ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji) 
     253            END DO 
     254         END DO 
     255 
     256         ! ghost cells: 
     257         ind1 = jpjglo - nbghostcells 
     258         ind2 = jpjglo 
     259         DO jj = mj0(ind1), mj1(ind2) 
     260            DO ji = 1, jpi 
     261               ztabramp(ji,jj) = zmsknorth(ji) 
     262            END DO 
     263         END DO 
    153264 
    154265      ENDIF 
    … …  
    156267      ! Tracers 
    157268      IF( .NOT. spongedoneT ) THEN 
    158          fsaht_spu(:,:) = 0._wp 
    159          fsaht_spv(:,:) = 0._wp 
    160          DO jj = 2, jpjm1 
    161             DO ji = 2, jpim1   ! vector opt. 
    162                fsaht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) 
    163                fsaht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) 
    164             END DO 
    165          END DO 
    166          CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spu, 'U', 1. )   ! Lateral boundary conditions 
    167          CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spv, 'V', 1. ) 
    168           
     269         fspu(:,:) = 0._wp 
     270         fspv(:,:) = 0._wp 
     271         DO_2D_00_00 
     272            fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) * ssumask(ji,jj) 
     273            fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) * ssvmask(ji,jj) 
     274         END_2D 
     275         CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. )   ! Lateral boundary conditions 
     276         CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. ) 
     277 
    169278         spongedoneT = .TRUE. 
    170279      ENDIF 
    … …  
    172281      ! Dynamics 
    173282      IF( .NOT. spongedoneU ) THEN 
    174          fsahm_spt(:,:) = 0._wp 
    175          fsahm_spf(:,:) = 0._wp 
    176          DO jj = 2, jpjm1 
    177             DO ji = 2, jpim1   ! vector opt. 
    178                fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj) 
    179                fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1) & 
    180                                                      &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) ) 
    181             END DO 
    182          END DO 
    183          CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spt, 'T', 1. )   ! Lateral boundary conditions 
    184          CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spf, 'F', 1. ) 
     283         fspt(:,:) = 0._wp 
     284         fspf(:,:) = 0._wp 
     285         DO_2D_00_00 
     286            fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 
     287            fspf(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1)   & 
     288                                  &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) ) & 
     289                                  &  * ssvmask(ji,jj) * ssvmask(ji,jj+1) 
     290         END_2D 
     291         CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. )   ! Lateral boundary conditions 
     292         CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. ) 
    185293          
    186294         spongedoneU = .TRUE. 
    187295      ENDIF 
     296 
     297#if defined key_vertical 
     298      ! Remove vertical interpolation where not needed: 
     299      DO_2D_00_00 
     300         IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 
     301         &   (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 
     302! 
     303         IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
     304         &   (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 
     305! 
     306         IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 
     307         &   (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 
     308! 
     309         IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 
     310         IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 
     311         IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 
     312      END_2D 
     313      ! 
     314      ztabramp(:,:) = REAL( mbkt_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. ) 
     315      mbkt_parent(:,:) = NINT( ztabramp(:,:) ) 
     316      ztabramp(:,:) = REAL( mbku_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. ) 
     317      mbku_parent(:,:) = NINT( ztabramp(:,:) ) 
     318      ztabramp(:,:) = REAL( mbkv_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. ) 
     319      mbkv_parent(:,:) = NINT( ztabramp(:,:) ) 
     320#endif 
    188321      ! 
    189322#endif 
    … …  
    201334      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    202335      INTEGER  ::   iku, ikv 
    203       REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
     336      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax 
    204337      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
    205338      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
    … …  
    210343      REAL(wp), DIMENSION(1:jpk) :: h_out 
    211344      INTEGER :: N_in, N_out 
    212       REAL(wp) :: h_diff 
    213345      !!---------------------------------------------------------------------- 
    214346      ! 
    … …  
    218350               DO jj=j1,j2 
    219351                  DO ji=i1,i2 
    220                      tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) 
     352                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a) 
    221353                  END DO 
    222354               END DO 
    … …  
    225357 
    226358# if defined key_vertical 
    227          DO jk=k1,k2 
    228             DO jj=j1,j2 
    229                DO ji=i1,i2 
    230                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
    231                END DO 
    232             END DO 
    233          END DO 
     359        ! Interpolate thicknesses 
     360        ! Warning: these are masked, hence extrapolated prior interpolation. 
     361        DO jk=k1,k2 
     362           DO jj=j1,j2 
     363              DO ji=i1,i2 
     364                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 
     365              END DO 
     366           END DO 
     367        END DO 
     368 
     369        ! Extrapolate thicknesses in partial bottom cells: 
     370        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     371        IF (ln_zps) THEN 
     372           DO jj=j1,j2 
     373              DO ji=i1,i2 
     374                  jk = mbkt(ji,jj) 
     375                  tabres(ji,jj,jk,jpts+1) = 0._wp 
     376              END DO 
     377           END DO            
     378        END IF 
     379      
     380        ! Save ssh at last level: 
     381        IF (.NOT.ln_linssh) THEN 
     382           tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)  
     383        ELSE 
     384           tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 
     385        END IF       
    234386# endif 
    235387 
    … …  
    237389         ! 
    238390# if defined key_vertical 
    239          tabres_child(:,:,:,:) = 0. 
     391 
     392         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 
     393 
    240394         DO jj=j1,j2 
    241395            DO ji=i1,i2 
    242                N_in = 0 
    243                DO jk=k1,k2 !k2 = jpk of parent grid 
    244                   IF (tabres(ji,jj,jk,n2) == 0) EXIT 
    245                   N_in = N_in + 1 
     396               tabres_child(ji,jj,:,:) = 0._wp  
     397               N_in = mbkt_parent(ji,jj) 
     398               zhtot = 0._wp 
     399               DO jk=1,N_in !k2 = jpk of parent grid 
     400                  IF (jk==N_in) THEN 
     401                     h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 
     402                  ELSE 
     403                     h_in(jk) = tabres(ji,jj,jk,n2) 
     404                  ENDIF 
     405                  zhtot = zhtot + h_in(jk) 
    246406                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    247                   h_in(N_in) = tabres(ji,jj,jk,n2) 
    248407               END DO 
    249408               N_out = 0 
    … …  
    251410                  IF (tmask(ji,jj,jk) == 0) EXIT  
    252411                  N_out = N_out + 1 
    253                   h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     412                  h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    254413               ENDDO 
    255                IF (N_in > 0) THEN 
    256                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    257                   tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 
    258                   DO jn=1,jpts 
    259                      call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
    260                   ENDDO 
     414 
     415               ! Account for small differences in free-surface 
     416               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     417                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     418               ELSE 
     419                  h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     420               ENDIF 
     421               IF (N_in*N_out > 0) THEN 
     422                  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) 
    261423               ENDIF 
    262424            ENDDO 
    … …  
    268430               DO jk=1,jpkm1 
    269431# if defined key_vertical 
    270                   tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts) 
     432                  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) 
    271433# else 
    272                   tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts) 
     434                  tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 
    273435# endif 
    274436               ENDDO 
    275437            ENDDO 
    276438         ENDDO 
     439 
     440         !* set relaxation time scale 
     441         IF( neuler == 0 .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_tra  / (        rdt ) 
     442         ELSE                                          ;   ztrelax =   rn_trelax_tra  / (2._wp * rdt ) 
     443         ENDIF 
    277444 
    278445         DO jn = 1, jpts             
    … …  
    281448               DO jj = j1,j2 
    282449                  DO ji = i1,i2-1 
    283                      zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 
     450                     zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    284451                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
    285452                  END DO 
    … …  
    288455               DO ji = i1,i2 
    289456                  DO jj = j1,j2-1 
    290                      zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 
     457                     zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 
    291458                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    292459                  END DO 
    … …  
    310477                  DO ji = i1+1,i2-1 
    311478                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN  
    312                         zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     479                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 
    313480                        ! horizontal diffusive trends 
    314                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) 
     481                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) & 
     482                             &  - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)  
    315483                        ! add it to the general tracer trends 
    316                         tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa 
     484                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa 
    317485                     ENDIF 
    318486                  END DO 
    … …  
    339507 
    340508      ! sponge parameters  
    341       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff 
     509      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 
    342510      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 
    343511      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
    … …  
    346514      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    347515      REAL(wp), DIMENSION(1:jpk) :: h_out 
    348       INTEGER ::N_in,N_out 
     516      INTEGER ::N_in, N_out 
    349517      !!---------------------------------------------     
    350518      ! 
    351519      IF( before ) THEN 
    352          DO jk=1,jpkm1 
     520         DO jk=k1,k2 
    353521            DO jj=j1,j2 
    354522               DO ji=i1,i2 
    355                   tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 
     523                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 
    356524# if defined key_vertical 
    357                   tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk) 
     525                  tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 
    358526# endif 
    359527               END DO 
    … …  
    361529         END DO 
    362530 
     531# if defined key_vertical 
     532         ! Extrapolate thicknesses in partial bottom cells: 
     533         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     534         IF (ln_zps) THEN 
     535            DO jj=j1,j2 
     536               DO ji=i1,i2 
     537                  jk = mbku(ji,jj) 
     538                  tabres(ji,jj,jk,m2) = 0._wp 
     539               END DO 
     540            END DO            
     541         END IF 
     542        ! Save ssh at last level: 
     543        tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     544        IF (.NOT.ln_linssh) THEN 
     545           ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 
     546           DO jk=1,jpk 
     547              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) 
     548           END DO 
     549           tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 
     550        END IF  
     551# endif 
     552 
    363553      ELSE 
    364554 
    365555# if defined key_vertical 
    366          tabres_child(:,:,:) = 0._wp 
     556         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     557 
    367558         DO jj=j1,j2 
    368559            DO ji=i1,i2 
    369                N_in = 0 
    370                DO jk=k1,k2 
    371                   IF (tabres(ji,jj,jk,m2) == 0) EXIT 
    372                   N_in = N_in + 1 
     560               tabres_child(ji,jj,:) = 0._wp 
     561               N_in = mbku_parent(ji,jj) 
     562               zhtot = 0._wp 
     563               DO jk=1,N_in 
     564                  IF (jk==N_in) THEN 
     565                     h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     566                  ELSE 
     567                     h_in(jk) = tabres(ji,jj,jk,m2) 
     568                  ENDIF 
     569                  zhtot = zhtot + h_in(jk) 
    373570                  tabin(jk) = tabres(ji,jj,jk,m1) 
    374                   h_in(N_in) = tabres(ji,jj,jk,m2) 
    375               ENDDO 
    376               ! 
    377               IF (N_in == 0) THEN 
    378                  tabres_child(ji,jj,:) = 0. 
    379                  CYCLE 
    380               ENDIF 
    381           
    382               N_out = 0 
    383               DO jk=1,jpk 
    384                  if (umask(ji,jj,jk) == 0) EXIT 
    385                  N_out = N_out + 1 
    386                  h_out(N_out) = e3u_n(ji,jj,jk) 
    387               ENDDO 
    388           
    389               IF (N_out == 0) THEN 
    390                  tabres_child(ji,jj,:) = 0. 
    391                  CYCLE 
    392               ENDIF 
    393           
    394               IF (N_in * N_out > 0) THEN 
    395                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    396                  if (h_diff < -1.e4) then 
    397                     print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
    398                  endif 
    399               ENDIF 
    400               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) 
    401           
     571               ENDDO 
     572               !          
     573               N_out = 0 
     574               DO jk=1,jpk 
     575                  IF (umask(ji,jj,jk) == 0) EXIT 
     576                  N_out = N_out + 1 
     577                  h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 
     578               ENDDO 
     579 
     580               ! Account for small differences in free-surface 
     581               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     582                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     583               ELSE 
     584                  h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     585               ENDIF 
     586                   
     587               IF (N_in * N_out > 0) THEN 
     588                  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) 
     589               ENDIF  
    402590            ENDDO 
    403591         ENDDO 
    404592 
    405          ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
     593         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
    406594#else 
    407          ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
     595         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
    408596#endif 
     597         !* set relaxation time scale 
     598         IF( neuler == 0 .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rdt ) 
     599         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rdt ) 
     600         ENDIF 
    409601         ! 
    410602         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    416608            DO jj = j1,j2 
    417609               DO ji = i1+1,i2   ! vector opt. 
    418                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
    419                   hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) & 
    420                                      &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 
     610                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj) 
     611                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) & 
     612                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr 
    421613               END DO 
    422614            END DO 
    … …  
    424616            DO jj = j1,j2-1 
    425617               DO ji = i1,i2   ! vector opt. 
    426                   zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj) 
     618                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 
    427619                  rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk)   & 
    428620                                    &   +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) ) * fmask(ji,jj,jk) * zbtr  
    … …  
    439631                     ze1v = hdivdiff(ji,jj,jk) 
    440632                     ! horizontal diffusive trends 
    441                      zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   & 
    442                            + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) 
     633                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   & 
     634                         & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &  
     635                         & - ztrelax  * fspu(ji,jj) * ubdiff(ji,jj,jk) 
    443636 
    444637                     ! add it to the general momentum trends 
    445                      ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    446  
     638                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua                                  
    447639                  END DO 
    448640               ENDIF 
    … …  
    465657 
    466658                     ! horizontal diffusive trends 
    467                      zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   & 
     659                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   & 
    468660                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj) 
    469661 
    470662                     ! add it to the general momentum trends 
    471                      va(ji,jj,jk) = va(ji,jj,jk) + zva 
     663                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva 
    472664                  END DO 
    473665               ENDIF 
    … …  
    492684      ! 
    493685      INTEGER  ::   ji, jj, jk, imax 
    494       REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff 
     686      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 
    495687      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 
    496688      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
    … …  
    503695       
    504696      IF( before ) THEN  
    505          DO jk=1,jpkm1 
     697         DO jk=k1,k2 
    506698            DO jj=j1,j2 
    507699               DO ji=i1,i2 
    508                   tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 
     700                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) 
    509701# if defined key_vertical 
    510                   tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk) 
     702                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a) 
    511703# endif 
    512704               END DO 
    513705            END DO 
    514706         END DO 
     707 
     708# if defined key_vertical 
     709         ! Extrapolate thicknesses in partial bottom cells: 
     710         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     711         IF (ln_zps) THEN 
     712            DO jj=j1,j2 
     713               DO ji=i1,i2 
     714                  jk = mbkv(ji,jj) 
     715                  tabres(ji,jj,jk,m2) = 0._wp 
     716               END DO 
     717            END DO            
     718         END IF 
     719        ! Save ssh at last level: 
     720        tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
     721        IF (.NOT.ln_linssh) THEN 
     722           ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 
     723           DO jk=1,jpk 
     724              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) 
     725           END DO 
     726           tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 
     727        END IF  
     728# endif 
     729 
    515730      ELSE 
    516731 
    517732# if defined key_vertical 
    518          tabres_child(:,:,:) = 0._wp 
     733         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 
    519734         DO jj=j1,j2 
    520735            DO ji=i1,i2 
    521                N_in = 0 
    522                DO jk=k1,k2 
    523                   IF (tabres(ji,jj,jk,m2) == 0) EXIT 
    524                   N_in = N_in + 1 
     736               tabres_child(ji,jj,:) = 0._wp 
     737               N_in = mbkv_parent(ji,jj) 
     738               zhtot = 0._wp 
     739               DO jk=1,N_in 
     740                  IF (jk==N_in) THEN 
     741                     h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 
     742                  ELSE 
     743                     h_in(jk) = tabres(ji,jj,jk,m2) 
     744                  ENDIF 
     745                  zhtot = zhtot + h_in(jk) 
    525746                  tabin(jk) = tabres(ji,jj,jk,m1) 
    526                   h_in(N_in) = tabres(ji,jj,jk,m2) 
    527               ENDDO 
     747               ENDDO 
     748               !           
     749               N_out = 0 
     750               DO jk=1,jpk 
     751                  IF (vmask(ji,jj,jk) == 0) EXIT 
     752                  N_out = N_out + 1 
     753                  h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 
     754               ENDDO 
     755 
     756               ! Account for small differences in free-surface 
     757               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 
     758                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 
     759               ELSE 
     760                  h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 
     761               ENDIF 
    528762          
    529               IF (N_in == 0) THEN 
    530                  tabres_child(ji,jj,:) = 0. 
    531                  CYCLE 
    532               ENDIF 
    533           
    534               N_out = 0 
    535               DO jk=1,jpk 
    536                  if (vmask(ji,jj,jk) == 0) EXIT 
    537                  N_out = N_out + 1 
    538                  h_out(N_out) = e3v_n(ji,jj,jk) 
    539               ENDDO 
    540           
    541               IF (N_in * N_out > 0) THEN 
    542                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    543                  if (h_diff < -1.e4) then 
    544                     print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
    545                  endif 
    546               ENDIF 
    547               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) 
     763               IF (N_in * N_out > 0) THEN 
     764                  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) 
     765               ENDIF 
    548766            ENDDO 
    549767         ENDDO 
    550768 
    551          vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
     769         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
    552770# else 
    553          vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
     771         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
    554772# endif 
     773         !* set relaxation time scale 
     774         IF( neuler == 0 .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rdt ) 
     775         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rdt ) 
     776         ENDIF 
    555777         ! 
    556778         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    562784            DO jj = j1+1,j2 
    563785               DO ji = i1,i2   ! vector opt. 
    564                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
    565                   hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  & 
    566                                      &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr 
     786                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj) 
     787                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  & 
     788                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr 
    567789               END DO 
    568790            END DO 
    569791            DO jj = j1,j2 
    570792               DO ji = i1,i2-1   ! vector opt. 
    571                   zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj) 
     793                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 
    572794                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &  
    573795                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr 
    … …  
    586808               IF( .NOT. tabspongedone_u(ji,jj) ) THEN 
    587809                  DO jk = 1, jpkm1 
    588                      ua(ji,jj,jk) = ua(ji,jj,jk)                                                               & 
    589                         & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  & 
     810                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                               & 
     811                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  & 
    590812                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj) 
    591813                  END DO 
    … …  
    600822               IF( .NOT. tabspongedone_v(ji,jj) ) THEN 
    601823                  DO jk = 1, jpkm1 
    602                      va(ji,jj,jk) = va(ji,jj,jk)                                                                  & 
    603                         &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   & 
    604                         &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj) 
     824                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                                  & 
     825                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   & 
     826                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                      & 
     827                        &  - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk) 
    605828                  END DO 
    606829               ENDIF 
Note: See TracChangeset for help on using the changeset viewer.