Ignore:
Timestamp:
2015-07-31T10:55:56+02:00 (6 years ago)
Author:
timgraham
Message:

Merge of AGRIF branch (branches/2014/dev_r4765_CNRS_agrif) onto the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r4153 r5656  
    11#define SPONGE && define SPONGE_TOP 
    22 
    3 Module agrif_opa_sponge 
     3MODULE agrif_opa_sponge 
    44#if defined key_agrif  && ! defined key_offline 
    55   USE par_oce 
     
    99   USE agrif_oce 
    1010   USE wrk_nemo   
     11   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1112 
    1213   IMPLICIT NONE 
    1314   PRIVATE 
    1415 
    15    PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn 
    16  
    17   !! * Substitutions 
     16   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
     17   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge 
     18 
     19   !! * Substitutions 
    1820#  include "domzgr_substitute.h90" 
    1921   !!---------------------------------------------------------------------- 
     
    2325   !!---------------------------------------------------------------------- 
    2426 
    25    CONTAINS 
     27CONTAINS 
    2628 
    2729   SUBROUTINE Agrif_Sponge_Tra 
     
    3032      !!--------------------------------------------- 
    3133      !! 
    32       INTEGER :: ji,jj,jk,jn 
    3334      REAL(wp) :: timecoeff 
    34       REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
    35       REAL(wp), POINTER, DIMENSION(:,:    ) :: ztu, ztv 
    36       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztab 
    37       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: tsbdiff 
    3835 
    3936#if defined SPONGE 
    40       CALL wrk_alloc( jpi, jpj, ztu, ztv ) 
    41       CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  ) 
    42  
    4337      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 
    4438 
     39      CALL Agrif_Sponge 
    4540      Agrif_SpecialValue=0. 
    4641      Agrif_UseSpecialValue = .TRUE. 
    47       ztab = 0.e0 
    48       CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn) 
     42      tabspongedone_tsn = .FALSE. 
     43 
     44      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge) 
     45 
    4946      Agrif_UseSpecialValue = .FALSE. 
    50  
    51       tsbdiff(:,:,:,:) = tsb(:,:,:,:) - ztab(:,:,:,:) 
    52  
    53       CALL Agrif_Sponge 
    54  
    55       DO jn = 1, jpts 
    56          DO jk = 1, jpkm1 
    57             ! 
    58             DO jj = 1, jpjm1 
    59                DO ji = 1, jpim1 
    60                   zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk) 
    61                   zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk) 
    62                   ztu(ji,jj) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    63                   ztv(ji,jj) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
    64                ENDDO 
    65             ENDDO 
    66  
    67             DO jj = 2, jpjm1 
    68                DO ji = 2, jpim1 
    69                   zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk) 
    70                   ! horizontal diffusive trends 
    71                   ztsa = zbtr * (  ztu(ji,jj) - ztu(ji-1,jj  )   & 
    72                   &              + ztv(ji,jj) - ztv(ji  ,jj-1)  ) 
    73                   ! add it to the general tracer trends 
    74                   tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa 
    75                END DO 
    76             END DO 
    77             ! 
    78          ENDDO 
    79       ENDDO 
    80  
    81       CALL wrk_dealloc( jpi, jpj, ztu, ztv ) 
    82       CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  ) 
    8347#endif 
    8448 
     
    9054      !!--------------------------------------------- 
    9155      !! 
    92       INTEGER :: ji,jj,jk 
    9356      REAL(wp) :: timecoeff 
    94       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
    95       REAL(wp), POINTER, DIMENSION(:,:,:) :: ubdiff, vbdiff 
    96       REAL(wp), POINTER, DIMENSION(:,:,:) :: rotdiff, hdivdiff 
    97       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab 
    9857 
    9958#if defined SPONGE 
    100       CALL wrk_alloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff ) 
    101  
    10259      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 
    10360 
    10461      Agrif_SpecialValue=0. 
    10562      Agrif_UseSpecialValue = ln_spc_dyn 
    106       ztab = 0.e0 
    107       CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun) 
     63 
     64      tabspongedone_u = .FALSE. 
     65      tabspongedone_v = .FALSE.          
     66      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge) 
     67 
     68      tabspongedone_u = .FALSE. 
     69      tabspongedone_v = .FALSE. 
     70      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge) 
     71 
    10872      Agrif_UseSpecialValue = .FALSE. 
    109  
    110       ubdiff(:,:,:) = ( ub(:,:,:) - ztab(:,:,:) ) * umask(:,:,:) 
    111  
    112       ztab = 0.e0 
    113       Agrif_SpecialValue=0. 
    114       Agrif_UseSpecialValue = ln_spc_dyn 
    115       CALL Agrif_Bc_Variable(ztab, va_id,calledweight=timecoeff,procname=interpvn) 
    116       Agrif_UseSpecialValue = .FALSE. 
    117  
    118       vbdiff(:,:,:) = ( vb(:,:,:) - ztab(:,:,:) ) * vmask(:,:,:) 
    119  
    120       CALL Agrif_Sponge 
    121  
    122       DO jk = 1,jpkm1 
    123          ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:) 
    124          vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:) 
    125       ENDDO 
    126        
    127       hdivdiff = 0. 
    128       rotdiff = 0. 
    129  
    130       DO jk = 1, jpkm1                                 ! Horizontal slab 
    131          !                                             ! =============== 
    132  
    133          !                                             ! -------- 
    134          ! Horizontal divergence                       !   div 
    135          !                                             ! -------- 
    136          DO jj = 2, jpjm1 
    137             DO ji = 2, jpim1   ! vector opt. 
    138                zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk) 
    139                hdivdiff(ji,jj,jk) =  (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * ubdiff(ji  ,jj  ,jk)     & 
    140                   &                   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * ubdiff(ji-1,jj  ,jk)     & 
    141                   &                   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vbdiff(ji  ,jj  ,jk)     & 
    142                   &                   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vbdiff(ji  ,jj-1,jk)  ) * zbtr 
    143             END DO 
    144          END DO 
    145  
    146          DO jj = 1, jpjm1 
    147             DO ji = 1, jpim1   ! vector opt. 
    148                zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk) 
    149                rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    & 
    150                   &                 - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) & 
    151                   &               * fmask(ji,jj,jk) * zbtr 
    152             END DO 
    153          END DO 
    154  
    155       ENDDO 
    156  
    157       !                                                ! =============== 
    158       DO jk = 1, jpkm1                                 ! Horizontal slab 
    159          !                                             ! =============== 
    160          DO jj = 2, jpjm1 
    161             DO ji = 2, jpim1   ! vector opt. 
    162                ! horizontal diffusive trends 
    163                zua = - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
    164                      + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk) ) / e1u(ji,jj) 
    165  
    166                zva = + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
    167                      + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) / e2v(ji,jj) 
    168                ! add it to the general momentum trends 
    169                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    170                va(ji,jj,jk) = va(ji,jj,jk) + zva 
    171             END DO 
    172          END DO 
    173          !                                             ! =============== 
    174       END DO                                           !   End of slab 
    175       !                                                ! =============== 
    176       CALL wrk_dealloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff ) 
    17773#endif 
    17874 
     
    19995         CALL wrk_alloc( jpi, jpj, ztabramp ) 
    20096 
    201          ispongearea  = 2 + 2 * Agrif_irhox() 
     97         ispongearea  = 2 + nn_sponge_len * Agrif_irhox() 
    20298         ilci = nlci - ispongearea 
    20399         ilcj = nlcj - ispongearea  
    204100         z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
    205          spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    206  
    207          ztabramp(:,:) = 0. 
     101 
     102         ztabramp(:,:) = 0._wp 
    208103 
    209104         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
     
    254149      ! Tracers 
    255150      IF( .NOT. spongedoneT ) THEN 
    256          spe1ur(:,:) = 0. 
    257          spe2vr(:,:) = 0. 
    258  
    259          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    260             spe1ur(2:ispongearea-1,:       ) = visc_tra                                        & 
    261                &                             *    0.5 * (  ztabramp(2:ispongearea-1,:      )   & 
    262                &                                         + ztabramp(3:ispongearea  ,:      ) ) & 
    263                &                             * e2u(2:ispongearea-1,:) / e1u(2:ispongearea-1,:) 
    264  
    265             spe2vr(2:ispongearea  ,1:jpjm1 ) = visc_tra                                        & 
    266                &                             *    0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1)   & 
    267                &                                         + ztabramp(2:ispongearea,2  :jpj  ) ) & 
    268                &                             * e1v(2:ispongearea,1:jpjm1) / e2v(2:ispongearea,1:jpjm1) 
    269          ENDIF 
    270  
    271          IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    272             spe1ur(ilci+1:nlci-2,:        ) = visc_tra                                   & 
    273                &                            * 0.5 * (  ztabramp(ilci+1:nlci-2,:      )   &  
    274                &                                     + ztabramp(ilci+2:nlci-1,:      ) ) & 
    275                &                            * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:) 
    276  
    277             spe2vr(ilci+1:nlci-1,1:jpjm1  )  = visc_tra                                  & 
    278                &                            * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1)   &  
    279                &                                     + ztabramp(ilci+1:nlci-1,2:jpj  ) ) &  
    280                &                            * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1) 
    281          ENDIF 
    282  
    283          IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    284             spe1ur(1:jpim1,2:ispongearea  ) = visc_tra                                     & 
    285                &                            * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  )   &  
    286                &                                     + ztabramp(2:jpi  ,2:ispongearea  ) ) & 
    287                &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea) 
    288     
    289             spe2vr(:      ,2:ispongearea-1) = visc_tra                                     & 
    290                &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   & 
    291                &                                     + ztabramp(:      ,3:ispongearea  ) ) & 
    292                &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1) 
    293          ENDIF 
    294  
    295          IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    296             spe1ur(1:jpim1,ilcj+1:nlcj-1) = visc_tra                                   & 
    297                &                          * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1)   & 
    298                &                                   + ztabramp(2:jpi  ,ilcj+1:nlcj-1) ) & 
    299                &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1) 
    300  
    301             spe2vr(:      ,ilcj+1:nlcj-2) = visc_tra                                   & 
    302                &                          * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2)   & 
    303                &                                   + ztabramp(:      ,ilcj+2:nlcj-1) ) & 
    304                &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2) 
    305          ENDIF 
     151         fsaht_spu(:,:) = 0._wp 
     152         fsaht_spv(:,:) = 0._wp 
     153         DO jj = 2, jpjm1 
     154            DO ji = 2, jpim1   ! vector opt. 
     155               fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj  )) 
     156               fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji  ,jj+1)) 
     157            END DO 
     158         END DO 
     159 
     160         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions 
     161         CALL lbc_lnk( fsaht_spv, 'V', 1. ) 
    306162         spongedoneT = .TRUE. 
    307163      ENDIF 
     
    309165      ! Dynamics 
    310166      IF( .NOT. spongedoneU ) THEN 
    311          spe1ur2(:,:) = 0. 
    312          spe2vr2(:,:) = 0. 
    313  
    314          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    315             spe1ur2(2:ispongearea-1,:      ) = visc_dyn                                   & 
    316                &                             * 0.5 * (  ztabramp(2:ispongearea-1,:      ) & 
    317                &                                      + ztabramp(3:ispongearea  ,:      ) ) 
    318             spe2vr2(2:ispongearea  ,1:jpjm1) = visc_dyn                                   & 
    319                &                             * 0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1) & 
    320                &                                      + ztabramp(2:ispongearea  ,2:jpj  ) )  
    321          ENDIF 
    322  
    323          IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    324             spe1ur2(ilci+1:nlci-2  ,:      ) = visc_dyn                                   & 
    325                &                             * 0.5 * (  ztabramp(ilci+1:nlci-2, :       ) & 
    326                &                                      + ztabramp(ilci+2:nlci-1, :       ) )                       
    327             spe2vr2(ilci+1:nlci-1  ,1:jpjm1) = visc_dyn                                   & 
    328                &                             * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1  ) & 
    329                &                                      + ztabramp(ilci+1:nlci-1,2:jpj    ) )  
    330          ENDIF 
    331  
    332          IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    333             spe1ur2(1:jpim1,2:ispongearea  ) = visc_dyn                                   &   
    334                &                             * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  ) & 
    335                &                                      + ztabramp(2:jpi  ,2:ispongearea  ) )  
    336             spe2vr2(:      ,2:ispongearea-1) = visc_dyn                                   & 
    337                &                             * 0.5 * (  ztabramp(:      ,2:ispongearea-1) & 
    338                &                                      + ztabramp(:      ,3:ispongearea  ) ) 
    339          ENDIF 
    340  
    341          IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    342             spe1ur2(1:jpim1,ilcj+1:nlcj-1  ) = visc_dyn                                   & 
    343                &                             * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1  ) & 
    344                &                                      + ztabramp(2:jpi  ,ilcj+1:nlcj-1  ) )  
    345             spe2vr2(:      ,ilcj+1:nlcj-2  ) = visc_dyn                                   & 
    346                &                             * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2  ) & 
    347                &                                      + ztabramp(:      ,ilcj+2:nlcj-1  ) ) 
    348          ENDIF 
     167         fsahm_spt(:,:) = 0._wp 
     168         fsahm_spf(:,:) = 0._wp 
     169         DO jj = 2, jpjm1 
     170            DO ji = 2, jpim1   ! vector opt. 
     171               fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj) 
     172               fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) & 
     173                                                     &  +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) 
     174            END DO 
     175         END DO 
     176 
     177         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions 
     178         CALL lbc_lnk( fsahm_spf, 'F', 1. ) 
    349179         spongedoneU = .TRUE. 
    350          spbtr3(:,:) = 1. / ( e1f(:,:) * e2f(:,:) ) 
    351180      ENDIF 
    352181      ! 
     
    357186   END SUBROUTINE Agrif_Sponge 
    358187 
    359    SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2) 
    360       !!--------------------------------------------- 
    361       !!   *** ROUTINE interptsn *** 
     188   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before) 
     189      !!--------------------------------------------- 
     190      !!   *** ROUTINE interptsn_sponge *** 
    362191      !!--------------------------------------------- 
    363192      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    364193      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    365  
    366       tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
    367  
    368    END SUBROUTINE interptsn 
    369  
    370    SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2) 
    371       !!--------------------------------------------- 
    372       !!   *** ROUTINE interpun *** 
    373       !!--------------------------------------------- 
     194      LOGICAL, INTENT(in) :: before 
     195 
     196 
     197      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     198      INTEGER  ::   iku, ikv 
     199      REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
     200      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv 
     201      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 
     202      ! 
     203      IF (before) THEN 
     204         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
     205      ELSE    
     206    
     207         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)     
     208         DO jn = 1, jpts             
     209            DO jk = 1, jpkm1 
     210               DO jj = j1,j2-1 
     211                  DO ji = i1,i2-1 
     212                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
     213                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
     214                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )  
     215                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     216                  ENDDO 
     217               ENDDO 
     218 
     219               IF( ln_zps ) THEN      ! set gradient at partial step level 
     220                  DO jj = j1,j2-1 
     221                     DO ji = i1,i2-1 
     222                        ! last level 
     223                        iku = mbku(ji,jj) 
     224                        ikv = mbkv(ji,jj) 
     225                        IF( iku == jk ) THEN 
     226                           ztu(ji,jj,jk) = 0._wp 
     227                        ENDIF 
     228                        IF( ikv == jk ) THEN 
     229                           ztv(ji,jj,jk) = 0._wp 
     230                        ENDIF 
     231                     END DO 
     232                  END DO 
     233               ENDIF 
     234            ENDDO 
     235 
     236            DO jk = 1, jpkm1 
     237               DO jj = j1+1,j2-1 
     238                  DO ji = i1+1,i2-1 
     239 
     240                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN  
     241                        zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk) 
     242                        ! horizontal diffusive trends 
     243                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) 
     244                        ! add it to the general tracer trends 
     245                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa 
     246                     ENDIF 
     247 
     248                  ENDDO 
     249               ENDDO 
     250 
     251            ENDDO 
     252         ENDDO 
     253 
     254         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     255 
     256      ENDIF 
     257 
     258   END SUBROUTINE interptsn_sponge 
     259 
     260   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before) 
     261      !!--------------------------------------------- 
     262      !!   *** ROUTINE interpun_sponge *** 
     263      !!---------------------------------------------     
    374264      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    375265      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    376  
    377       tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2) 
    378  
    379    END SUBROUTINE interpun 
    380  
    381    SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2) 
    382       !!--------------------------------------------- 
    383       !!   *** ROUTINE interpvn *** 
    384       !!--------------------------------------------- 
     266      LOGICAL, INTENT(in) :: before 
     267 
     268      INTEGER :: ji,jj,jk 
     269 
     270      ! sponge parameters  
     271      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     272      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 
     273      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     274      INTEGER :: jmax 
     275      ! 
     276 
     277 
     278      IF (before) THEN 
     279         tabres = un(i1:i2,j1:j2,:) 
     280      ELSE 
     281 
     282         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 
     283 
     284         DO jk = 1, jpkm1                                 ! Horizontal slab 
     285            !                                             ! =============== 
     286 
     287            !                                             ! -------- 
     288            ! Horizontal divergence                       !   div 
     289            !                                             ! -------- 
     290            DO jj = j1,j2 
     291               DO ji = i1+1,i2   ! vector opt. 
     292                  zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
     293                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*fse3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) & 
     294                                     &   -e2u(ji-1,jj)*fse3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 
     295               END DO 
     296            END DO 
     297 
     298            DO jj = j1,j2-1 
     299               DO ji = i1,i2   ! vector opt. 
     300                  zbtr = r1_e12f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj) 
     301                  rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) & 
     302                                       +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) &  
     303                                    & ) * fmask(ji,jj,jk) * zbtr  
     304               END DO 
     305            END DO 
     306         ENDDO 
     307 
     308         ! 
     309 
     310 
     311 
     312         DO jj = j1+1, j2-1 
     313            DO ji = i1+1, i2-1   ! vector opt. 
     314 
     315               IF (.NOT. tabspongedone_u(ji,jj)) THEN 
     316                  DO jk = 1, jpkm1                                 ! Horizontal slab 
     317                     ze2u = rotdiff (ji,jj,jk) 
     318                     ze1v = hdivdiff(ji,jj,jk) 
     319                     ! horizontal diffusive trends 
     320                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u_n(ji,jj,jk) )   & 
     321                           + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / e1u(ji,jj) 
     322 
     323                     ! add it to the general momentum trends 
     324                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
     325 
     326                  END DO 
     327               ENDIF 
     328 
     329            END DO 
     330         END DO 
     331 
     332         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     333 
     334         jmax = j2-1 
     335         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3) 
     336 
     337         DO jj = j1+1, jmax 
     338            DO ji = i1+1, i2   ! vector opt. 
     339 
     340               IF (.NOT. tabspongedone_v(ji,jj)) THEN 
     341                  DO jk = 1, jpkm1                                 ! Horizontal slab 
     342                     ze2u = rotdiff (ji,jj,jk) 
     343                     ze1v = hdivdiff(ji,jj,jk) 
     344 
     345                     ! horizontal diffusive trends 
     346                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v_n(ji,jj,jk) )   & 
     347                           + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / e2v(ji,jj) 
     348 
     349                     ! add it to the general momentum trends 
     350                     va(ji,jj,jk) = va(ji,jj,jk) + zva 
     351                  END DO 
     352               ENDIF 
     353 
     354            END DO 
     355         END DO 
     356 
     357 
     358         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE. 
     359 
     360      ENDIF 
     361 
     362 
     363   END SUBROUTINE interpun_sponge 
     364 
     365 
     366   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir) 
     367      !!--------------------------------------------- 
     368      !!   *** ROUTINE interpvn_sponge *** 
     369      !!---------------------------------------------  
    385370      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    386371      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    387  
    388       tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2) 
    389  
    390    END SUBROUTINE interpvn 
     372      LOGICAL, INTENT(in) :: before 
     373      INTEGER, INTENT(in) :: nb , ndir 
     374 
     375      INTEGER :: ji,jj,jk 
     376 
     377      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     378 
     379      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff 
     380      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     381      INTEGER :: imax 
     382      ! 
     383 
     384      IF (before) THEN  
     385         tabres = vn(i1:i2,j1:j2,:) 
     386      ELSE 
     387 
     388         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:) 
     389 
     390         DO jk = 1, jpkm1                                 ! Horizontal slab 
     391            !                                             ! =============== 
     392 
     393            !                                             ! -------- 
     394            ! Horizontal divergence                       !   div 
     395            !                                             ! -------- 
     396            DO jj = j1+1,j2 
     397               DO ji = i1,i2   ! vector opt. 
     398                  zbtr = r1_e12t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
     399                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  & 
     400                                     &  -e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr 
     401               END DO 
     402            END DO 
     403            DO jj = j1,j2 
     404               DO ji = i1,i2-1   ! vector opt. 
     405                  zbtr = r1_e12f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj) 
     406                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &  
     407                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk) & 
     408                                    & ) * fmask(ji,jj,jk) * zbtr 
     409               END DO 
     410            END DO 
     411         ENDDO 
     412 
     413         !                                                ! =============== 
     414         !                                                 
     415 
     416         imax = i2-1 
     417         IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3) 
     418 
     419         DO jj = j1+1, j2 
     420            DO ji = i1+1, imax   ! vector opt. 
     421               IF (.NOT. tabspongedone_u(ji,jj)) THEN 
     422                  DO jk = 1, jpkm1                                 ! Horizontal slab 
     423                     ze2u = rotdiff (ji,jj,jk) 
     424                     ze1v = hdivdiff(ji,jj,jk) 
     425                     ! horizontal diffusive trends 
     426                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u_n(ji,jj,jk) ) + ( hdivdiff(ji+1,jj,jk) - ze1v) & 
     427                           / e1u(ji,jj) 
     428 
     429 
     430                     ! add it to the general momentum trends 
     431                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
     432                  END DO 
     433 
     434               ENDIF 
     435            END DO 
     436         END DO 
     437 
     438         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE. 
     439 
     440         DO jj = j1+1, j2-1 
     441            DO ji = i1+1, i2-1   ! vector opt. 
     442               IF (.NOT. tabspongedone_v(ji,jj)) THEN 
     443                  DO jk = 1, jpkm1                                 ! Horizontal slab 
     444                     ze2u = rotdiff (ji,jj,jk) 
     445                     ze1v = hdivdiff(ji,jj,jk) 
     446                     ! horizontal diffusive trends 
     447 
     448                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v_n(ji,jj,jk) ) + ( hdivdiff(ji,jj+1,jk) - ze1v) & 
     449                           / e2v(ji,jj) 
     450 
     451                     ! add it to the general momentum trends 
     452                     va(ji,jj,jk) = va(ji,jj,jk) + zva 
     453                  END DO 
     454               ENDIF 
     455            END DO 
     456         END DO 
     457         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     458      ENDIF 
     459 
     460   END SUBROUTINE interpvn_sponge 
    391461 
    392462#else 
Note: See TracChangeset for help on using the changeset viewer.