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 4785 for branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 – NEMO

Ignore:
Timestamp:
2014-09-24T14:03:02+02:00 (10 years ago)
Author:
rblod
Message:

dev_r4765_CNRS_agrif: First update of AGRIF for dynamic only (_flt and _ts), see ticket #1380 and associated wiki page

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r4153 r4785  
    1313   PRIVATE 
    1414 
    15    PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn 
     15   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn 
     16   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge 
    1617 
    1718  !! * Substitutions 
    … …  
    3839 
    3940#if defined SPONGE 
    40       CALL wrk_alloc( jpi, jpj, ztu, ztv ) 
    41       CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  ) 
    42  
    4341      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 
    44  
     42       
     43      CALL Agrif_Sponge 
    4544      Agrif_SpecialValue=0. 
    4645      Agrif_UseSpecialValue = .TRUE. 
    47       ztab = 0.e0 
    48       CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn) 
     46      tabspongedone = .FALSE. 
     47 
     48      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge) 
     49 
    4950      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  ) 
    8351#endif 
    8452 
    … …  
    9058      !!--------------------------------------------- 
    9159      !! 
    92       INTEGER :: ji,jj,jk 
    9360      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 
    9861 
    9962#if defined SPONGE 
    100       CALL wrk_alloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff ) 
    101  
    10263      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 
    10364 
    104       Agrif_SpecialValue=0. 
    105       Agrif_UseSpecialValue = ln_spc_dyn 
    106       ztab = 0.e0 
    107       CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun) 
    108       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 ) 
     65         Agrif_SpecialValue=0. 
     66         Agrif_UseSpecialValue = ln_spc_dyn 
     67 
     68         tabspongedone_u = .FALSE. 
     69         tabspongedone_v = .FALSE.          
     70         CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge) 
     71 
     72         tabspongedone_u = .FALSE. 
     73         tabspongedone_v = .FALSE. 
     74         CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge) 
     75 
     76         Agrif_UseSpecialValue = .FALSE. 
    17777#endif 
    17878 
    … …  
    18585      INTEGER  :: ji,jj,jk 
    18686      INTEGER  :: ispongearea, ilci, ilcj 
    187       LOGICAL  :: ll_spdone 
    188       REAL(wp) :: z1spongearea, zramp 
    189       REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp 
     87      REAL(wp) :: z1spongearea 
     88      REAL(wp), POINTER, DIMENSION(:,:) :: zlocalviscsponge 
    19089 
    19190#if defined SPONGE || defined SPONGE_TOP 
    192       ll_spdone=.TRUE. 
    193       IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
    194          ! Define ramp from boundaries towards domain interior 
    195          ! at T-points 
    196          ! Store it in ztabramp 
    197          ll_spdone=.FALSE. 
    198  
    199          CALL wrk_alloc( jpi, jpj, ztabramp ) 
    200  
    201          ispongearea  = 2 + 2 * Agrif_irhox() 
    202          ilci = nlci - ispongearea 
    203          ilcj = nlcj - ispongearea  
    204          z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
    205          spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    206  
    207          ztabramp(:,:) = 0. 
    208  
    209          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    210             DO jj = 1, jpj 
    211                IF ( umask(2,jj,1) == 1._wp ) THEN 
    212                  DO ji = 2, ispongearea                   
    213                     ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea 
    214                  END DO 
    215                ENDIF 
    216             ENDDO 
    217          ENDIF 
    218  
    219          IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    220             DO jj = 1, jpj 
    221                IF ( umask(nlci-2,jj,1) == 1._wp ) THEN 
    222                   DO ji = ilci+1,nlci-1 
    223                      zramp = (ji - (ilci+1) ) * z1spongearea 
    224                      ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
    225                   ENDDO 
    226                ENDIF 
    227             ENDDO 
    228          ENDIF 
    229  
    230          IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    231             DO ji = 1, jpi 
    232                IF ( vmask(ji,2,1) == 1._wp ) THEN 
    233                   DO jj = 2, ispongearea 
    234                      zramp = ( ispongearea-jj ) * z1spongearea 
    235                      ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
    236                   END DO 
    237                ENDIF 
    238             ENDDO 
    239          ENDIF 
    240  
    241          IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    242             DO ji = 1, jpi 
    243                IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN 
    244                   DO jj = ilcj+1,nlcj-1 
    245                      zramp = (jj - (ilcj+1) ) * z1spongearea 
    246                      ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
    247                   END DO 
    248                ENDIF 
    249             ENDDO 
    250          ENDIF 
    251  
    252       ENDIF 
     91 
     92      CALL wrk_alloc( jpi, jpj, zlocalviscsponge ) 
     93 
     94      ispongearea  = 2 + 2 * Agrif_irhox() 
     95      ilci = nlci - ispongearea 
     96      ilcj = nlcj - ispongearea  
     97      z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
     98      spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    25399 
    254100      ! Tracers 
    255101      IF( .NOT. spongedoneT ) THEN 
     102         zlocalviscsponge(:,:) = 0. 
    256103         spe1ur(:,:) = 0. 
    257104         spe2vr(:,:) = 0. 
    258105 
    259106         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) 
     107            DO ji = 2, ispongearea 
     108               zlocalviscsponge(ji,:) = visc_tra * ( ispongearea-ji ) * z1spongearea 
     109            ENDDO 
     110            spe1ur(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      )   & 
     111               &                         +            zlocalviscsponge(3:ispongearea  ,:      ) ) & 
     112               &                         * e2u(2:ispongearea-1,:      ) / e1u(2:ispongearea-1,:      ) 
     113            spe2vr(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1)   & 
     114               &                         +            zlocalviscsponge(2:ispongearea,2  :jpj  ) ) & 
     115               &                         * e1v(2:ispongearea  ,1:jpjm1) / e2v(2:ispongearea  ,1:jpjm1) 
    269116         ENDIF 
    270117 
    271118         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) 
     119            DO ji = ilci+1,nlci-1 
     120               zlocalviscsponge(ji,:) = visc_tra * (ji - (ilci+1) ) * z1spongearea 
     121            ENDDO 
     122   
     123            spe1ur(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:)    &  
     124               &                          +          zlocalviscsponge(ilci+2:nlci-1,:) )  & 
     125               &                          * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:) 
     126 
     127            spe2vr(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1)    &  
     128               &                            +        zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) &  
     129               &                                   * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1) 
    281130         ENDIF 
    282131 
    283132         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  ) ) & 
     133            DO jj = 2, ispongearea 
     134               zlocalviscsponge(:,jj) = visc_tra * ( ispongearea-jj ) * z1spongearea 
     135            ENDDO 
     136            spe1ur(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea  ) &  
     137               &                            +         zlocalviscsponge(2:jpi  ,2:ispongearea) ) & 
    287138               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea) 
    288139    
    289             spe2vr(:      ,2:ispongearea-1) = visc_tra                                     & 
    290                &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   & 
    291                &                                     + ztabramp(:      ,3:ispongearea  ) ) & 
     140            spe2vr(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)       & 
     141               &                            +         zlocalviscsponge(:,3:ispongearea  )     ) & 
    292142               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1) 
    293143         ENDIF 
    294144 
    295145         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) ) & 
     146            DO jj = ilcj+1,nlcj-1 
     147               zlocalviscsponge(:,jj) = visc_tra * (jj - (ilcj+1) ) * z1spongearea 
     148            ENDDO 
     149            spe1ur(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1)   & 
     150               &                          +         zlocalviscsponge(2:jpi  ,ilcj+1:nlcj-1) ) & 
    299151               &                                * 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) ) & 
     152            spe2vr(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      )   & 
     153               &                          +         zlocalviscsponge(:,ilcj+2:nlcj-1)     )   & 
    304154               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2) 
    305155         ENDIF 
    … …  
    309159      ! Dynamics 
    310160      IF( .NOT. spongedoneU ) THEN 
     161         zlocalviscsponge(:,:) = 0. 
    311162         spe1ur2(:,:) = 0. 
    312163         spe2vr2(:,:) = 0. 
    313164 
    314165         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  ) )  
     166            DO ji = 2, ispongearea 
     167               zlocalviscsponge(ji,:) = visc_dyn * ( ispongearea-ji ) * z1spongearea 
     168            ENDDO 
     169            spe1ur2(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) & 
     170                                             &     +   zlocalviscsponge(3:ispongearea,:    ) ) 
     171            spe2vr2(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) & 
     172                                             &     +   zlocalviscsponge(2:ispongearea,2:jpj) )  
    321173         ENDIF 
    322174 
    323175         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    ) )  
     176            DO ji = ilci+1,nlci-1 
     177               zlocalviscsponge(ji,:) = visc_dyn * (ji - (ilci+1) ) * z1spongearea 
     178            ENDDO 
     179            spe1ur2(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) & 
     180                                           &        + zlocalviscsponge(ilci+2:nlci-1,:) )   
     181            spe2vr2(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) & 
     182                                           &        + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  )  
    330183         ENDIF 
    331184 
    332185         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  ) ) 
     186            DO jj = 2, ispongearea 
     187               zlocalviscsponge(:,jj) = visc_dyn * ( ispongearea-jj ) * z1spongearea 
     188            ENDDO 
     189            spe1ur2(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) & 
     190                                             &      + zlocalviscsponge(2:jpi,2:ispongearea) )  
     191            spe2vr2(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     & 
     192                                             &      + zlocalviscsponge(:,3:ispongearea)     ) 
    339193         ENDIF 
    340194 
    341195         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  ) ) 
     196            DO jj = ilcj+1,nlcj-1 
     197               zlocalviscsponge(:,jj) = visc_dyn * (jj - (ilcj+1) ) * z1spongearea 
     198            ENDDO 
     199            spe1ur2(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) & 
     200                                           &         + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) )  
     201            spe2vr2(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) & 
     202                                           &         + zlocalviscsponge(:,ilcj+2:nlcj-1)     ) 
    348203         ENDIF 
    349204         spongedoneU = .TRUE. 
    … …  
    351206      ENDIF 
    352207      ! 
    353       IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp ) 
     208      CALL wrk_dealloc( jpi, jpj, zlocalviscsponge ) 
    354209      ! 
    355210#endif 
    … …  
    357212   END SUBROUTINE Agrif_Sponge 
    358213 
    359    SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2) 
    360       !!--------------------------------------------- 
    361       !!   *** ROUTINE interptsn *** 
     214   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before) 
     215      !!--------------------------------------------- 
     216      !!   *** ROUTINE interptsn_sponge *** 
    362217      !!--------------------------------------------- 
    363218      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    364219      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       !!--------------------------------------------- 
     220      LOGICAL, INTENT(in) :: before 
     221        
     222      
     223      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     224 
     225      REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
     226      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ztu, ztv 
     227      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 
     228      ! 
     229          
     230          
     231         IF (before) THEN 
     232            tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
     233         ELSE       
     234 
     235            tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)       
     236            DO jn = 1, jpts 
     237               DO jk = 1, jpkm1 
     238                  
     239                  DO jj = j1,j2-1 
     240                     DO ji = i1,i2-1 
     241                        zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk) 
     242                        zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk) 
     243                        ztu(ji,jj) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     244                        ztv(ji,jj) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
     245                     ENDDO 
     246                  ENDDO 
     247                   
     248                  DO jj = j1+1,j2-1 
     249                     DO ji = i1+1,i2-1 
     250                         
     251                        if (.not. tabspongedone(ji,jj)) then  
     252                           zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk) 
     253            ! horizontal diffusive trends 
     254                           ztsa = zbtr * (  ztu(ji,jj) - ztu(ji-1,jj  ) + ztv(ji,jj) - ztv(ji  ,jj-1)  ) 
     255            ! add it to the general tracer trends 
     256                           tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa 
     257                         endif  
     258  
     259                       ENDDO 
     260                    ENDDO 
     261                     
     262                ENDDO 
     263             ENDDO 
     264              
     265             tabspongedone(i1+1:i2-1,j1+1:j2-1) = .TRUE. 
     266                          
     267    ENDIF 
     268                 
     269   END SUBROUTINE interptsn_sponge 
     270 
     271   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before) 
     272      !!--------------------------------------------- 
     273      !!   *** ROUTINE interpun_sponge *** 
     274      !!---------------------------------------------     
    374275      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    375276      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       !!--------------------------------------------- 
     277      LOGICAL, INTENT(in) :: before 
     278 
     279      INTEGER :: ji,jj,jk 
     280 
     281   ! sponge parameters  
     282      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     283      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 
     284      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     285      INTEGER :: jmax 
     286   ! 
     287       
     288 
     289      IF (before) THEN 
     290           
     291          tabres = un(i1:i2,j1:j2,:) 
     292 
     293      ELSE 
     294          
     295         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 
     296          
     297         DO jk=1,jpkm1 
     298            ubdiff(i1:i2,j1:j2,jk) = ubdiff(i1:i2,j1:j2,jk) * spe1ur2(i1:i2,j1:j2) 
     299         ENDDO 
     300 
     301         DO jk = 1, jpkm1                                 ! Horizontal slab 
     302!                                             ! =============== 
     303 
     304!                                             ! -------- 
     305! Horizontal divergence                       !   div 
     306!                                             ! -------- 
     307            DO jj = j1,j2 
     308               DO ji = i1+1,i2   ! vector opt. 
     309                  zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk) 
     310                  hdivdiff(ji,jj,jk) = (e2u(ji,jj)*fse3u(ji,jj,jk) * ubdiff(ji,jj,jk) - e2u(ji-1,jj)* fse3u(ji-1,jj  ,jk)  & 
     311                                       * ubdiff(ji-1,jj  ,jk) ) * zbtr 
     312               END DO 
     313            END DO 
     314 
     315            DO jj = j1,j2-1 
     316               DO ji = i1,i2   ! vector opt. 
     317                  zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk) 
     318                  rotdiff(ji,jj,jk) = (- e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) & 
     319                                      * fmask(ji,jj,jk) * zbtr  
     320               END DO 
     321            END DO 
     322         ENDDO 
     323 
     324! 
     325 
     326 
     327 
     328            DO jj = j1+1, j2-1 
     329               DO ji = i1+1, i2-1   ! vector opt. 
     330                   
     331                  if (.not. tabspongedone_u(ji,jj)) then 
     332                     DO jk = 1, jpkm1                                 ! Horizontal slab 
     333                        ze2u = rotdiff (ji,jj,jk) 
     334                        ze1v = hdivdiff(ji,jj,jk) 
     335! horizontal diffusive trends 
     336                        zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
     337                        + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / e1u(ji,jj) 
     338 
     339! add it to the general momentum trends 
     340                        ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
     341 
     342                     END DO                   
     343                  endif  
     344 
     345               END DO             
     346            END DO 
     347                   
     348            tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .true. 
     349  
     350         jmax = j2-1 
     351         If ((nbondj == 1).OR.(nbondj == 2)) jmax = min(jmax,nlcj-3) 
     352                                 
     353            DO jj = j1+1, jmax 
     354               DO ji = i1+1, i2   ! vector opt. 
     355                   
     356                  if (.not. tabspongedone_v(ji,jj)) then 
     357                     DO jk = 1, jpkm1                                 ! Horizontal slab 
     358                        ze2u = rotdiff (ji,jj,jk) 
     359                        ze1v = hdivdiff(ji,jj,jk) 
     360                         
     361! horizontal diffusive trends 
     362                        zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
     363                        + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / e2v(ji,jj) 
     364 
     365! add it to the general momentum trends 
     366                        va(ji,jj,jk) = va(ji,jj,jk) + zva 
     367                     END DO                   
     368                  endif  
     369 
     370               END DO             
     371            END DO 
     372 
     373             
     374            tabspongedone_v(i1+1:i2,j1+1:jmax) = .true. 
     375             
     376      ENDIF 
     377          
     378           
     379   END SUBROUTINE interpun_sponge 
     380  
     381    
     382   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir) 
     383     !!--------------------------------------------- 
     384      !!   *** ROUTINE interpvn_sponge *** 
     385      !!---------------------------------------------  
    385386      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    386387      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 
     388      LOGICAL, INTENT(in) :: before 
     389      INTEGER, INTENT(in) :: nb , ndir 
     390 
     391      INTEGER :: ji,jj,jk 
     392 
     393      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
     394       
     395      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff 
     396      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     397      INTEGER :: imax 
     398   ! 
     399      
     400      IF (before) THEN  
     401        tabres = vn(i1:i2,j1:j2,:) 
     402      ELSE 
     403 
     404         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:) 
     405             
     406         DO jk=1,jpkm1 
     407            vbdiff(i1:i2,j1:j2,jk) = vbdiff(i1:i2,j1:j2,jk) * spe2vr2(i1:i2,j1:j2) 
     408         ENDDO 
     409 
     410         DO jk = 1, jpkm1                                 ! Horizontal slab 
     411!                                             ! =============== 
     412 
     413!                                             ! -------- 
     414! Horizontal divergence                       !   div 
     415!                                             ! -------- 
     416            DO jj = j1+1,j2 
     417               DO ji = i1,i2   ! vector opt. 
     418                  zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk) 
     419                  hdivdiff(ji,jj,jk) = (e1v(ji,jj) * fse3v(ji,jj,jk) * vbdiff(ji,jj,jk) - e1v(ji  ,jj-1) & 
     420                                       * fse3v(ji  ,jj-1,jk)  * vbdiff(ji  ,jj-1,jk)  ) * zbtr 
     421               END DO 
     422            END DO 
     423 
     424            DO jj = j1,j2 
     425               DO ji = i1,i2-1   ! vector opt. 
     426                  zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk) 
     427                  rotdiff(ji,jj,jk) = (e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)) & 
     428                                      * fmask(ji,jj,jk) * zbtr 
     429               END DO 
     430            END DO 
     431 
     432         ENDDO 
     433 
     434!                                                ! =============== 
     435!                                                 
     436          
     437         imax = i2-1 
     438         If ((nbondi == 1).OR.(nbondi == 2)) imax = min(imax,nlci-3) 
     439                             
     440            DO jj = j1+1, j2 
     441               DO ji = i1+1, imax   ! vector opt. 
     442                  if (.not. tabspongedone_u(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                        zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) + ( hdivdiff(ji+1,jj,jk) - ze1v) & 
     448                        / e1u(ji,jj) 
     449 
     450 
     451! add it to the general momentum trends 
     452                       ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
     453                     END DO 
     454 
     455                   endif 
     456            END DO             
     457         END DO   
     458   
     459         tabspongedone_u(i1+1:imax,j1+1:j2) = .true. 
     460          
     461            DO jj = j1+1, j2-1 
     462               DO ji = i1+1, i2-1   ! vector opt. 
     463                  if (.not. tabspongedone_v(ji,jj)) then 
     464                     DO jk = 1, jpkm1                                 ! Horizontal slab 
     465                        ze2u = rotdiff (ji,jj,jk) 
     466                        ze1v = hdivdiff(ji,jj,jk) 
     467! horizontal diffusive trends 
     468 
     469                        zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) + ( hdivdiff(ji,jj+1,jk) - ze1v) & 
     470                        / e2v(ji,jj) 
     471 
     472! add it to the general momentum trends 
     473                       va(ji,jj,jk) = va(ji,jj,jk) + zva 
     474                     END DO 
     475 
     476                   endif 
     477            END DO             
     478         END DO           
     479          
     480         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .true. 
     481          
     482      ENDIF 
     483      
     484   END SUBROUTINE interpvn_sponge 
    391485 
    392486#else 
Note: See TracChangeset for help on using the changeset viewer.