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 4009 – NEMO

Changeset 4009


Ignore:
Timestamp:
2013-08-28T16:18:17+02:00 (10 years ago)
Author:
cbricaud
Message:

Correct Agrif sponge coefficients near bdy corners

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r3918 r4009  
    185185      INTEGER  :: ji,jj,jk 
    186186      INTEGER  :: ispongearea, ilci, ilcj 
    187       REAL(wp) :: z1spongearea 
    188       REAL(wp), POINTER, DIMENSION(:,:) :: zlocalviscsponge 
     187      LOGICAL  :: ll_spdone 
     188      REAL(wp) :: z1spongearea, zramp 
     189      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp 
    189190 
    190191#if defined SPONGE || defined SPONGE_TOP 
    191  
    192       CALL wrk_alloc( jpi, jpj, zlocalviscsponge ) 
    193  
    194       ispongearea  = 2 + 2 * Agrif_irhox() 
    195       ilci = nlci - ispongearea 
    196       ilcj = nlcj - ispongearea  
    197       z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
    198       spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     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 
    199253 
    200254      ! Tracers 
    201255      IF( .NOT. spongedoneT ) THEN 
    202          zlocalviscsponge(:,:) = 0. 
    203256         spe1ur(:,:) = 0. 
    204257         spe2vr(:,:) = 0. 
    205258 
    206259         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    207             DO ji = 2, ispongearea 
    208                zlocalviscsponge(ji,:) = visc_tra * ( ispongearea-ji ) * z1spongearea 
    209             ENDDO 
    210             spe1ur(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      )   & 
    211                &                         +            zlocalviscsponge(3:ispongearea  ,:      ) ) & 
    212                &                         * e2u(2:ispongearea-1,:      ) / e1u(2:ispongearea-1,:      ) 
    213             spe2vr(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1)   & 
    214                &                         +            zlocalviscsponge(2:ispongearea,2  :jpj  ) ) & 
    215                &                         * e1v(2:ispongearea  ,1:jpjm1) / e2v(2:ispongearea  ,1:jpjm1) 
     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) 
    216269         ENDIF 
    217270 
    218271         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    219             DO ji = ilci+1,nlci-1 
    220                zlocalviscsponge(ji,:) = visc_tra * (ji - (ilci+1) ) * z1spongearea 
    221             ENDDO 
    222    
    223             spe1ur(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:)    &  
    224                &                          +          zlocalviscsponge(ilci+2:nlci-1,:) )  & 
    225                &                          * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:) 
    226  
    227             spe2vr(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1)    &  
    228                &                            +        zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) &  
    229                &                                   * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1) 
     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) 
    230281         ENDIF 
    231282 
    232283         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    233             DO jj = 2, ispongearea 
    234                zlocalviscsponge(:,jj) = visc_tra * ( ispongearea-jj ) * z1spongearea 
    235             ENDDO 
    236             spe1ur(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea  ) &  
    237                &                            +         zlocalviscsponge(2:jpi  ,2:ispongearea) ) & 
     284            spe1ur(1:jpim1,2:ispongearea  ) = visc_tra                                     & 
     285               &                            * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  )   &  
     286               &                                     + ztabramp(2:jpi  ,2:ispongearea  ) ) & 
    238287               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea) 
    239288    
    240             spe2vr(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)       & 
    241                &                            +         zlocalviscsponge(:,3:ispongearea  )     ) & 
     289            spe2vr(:      ,2:ispongearea-1) = visc_tra                                     & 
     290               &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   & 
     291               &                                     + ztabramp(:      ,3:ispongearea  ) ) & 
    242292               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1) 
    243293         ENDIF 
    244294 
    245295         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    246             DO jj = ilcj+1,nlcj-1 
    247                zlocalviscsponge(:,jj) = visc_tra * (jj - (ilcj+1) ) * z1spongearea 
    248             ENDDO 
    249             spe1ur(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1)   & 
    250                &                          +         zlocalviscsponge(2:jpi  ,ilcj+1:nlcj-1) ) & 
     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) ) & 
    251299               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1) 
    252             spe2vr(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      )   & 
    253                &                          +         zlocalviscsponge(:,ilcj+2: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) ) & 
    254304               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2) 
    255305         ENDIF 
    … …  
    259309      ! Dynamics 
    260310      IF( .NOT. spongedoneU ) THEN 
    261          zlocalviscsponge(:,:) = 0. 
    262311         spe1ur2(:,:) = 0. 
    263312         spe2vr2(:,:) = 0. 
    264313 
    265314         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    266             DO ji = 2, ispongearea 
    267                zlocalviscsponge(ji,:) = visc_dyn * ( ispongearea-ji ) * z1spongearea 
    268             ENDDO 
    269             spe1ur2(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) & 
    270                                              &     +   zlocalviscsponge(3:ispongearea,:    ) ) 
    271             spe2vr2(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) & 
    272                                              &     +   zlocalviscsponge(2:ispongearea,2:jpj) )  
     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  ) )  
    273321         ENDIF 
    274322 
    275323         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    276             DO ji = ilci+1,nlci-1 
    277                zlocalviscsponge(ji,:) = visc_dyn * (ji - (ilci+1) ) * z1spongearea 
    278             ENDDO 
    279             spe1ur2(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) & 
    280                                            &        + zlocalviscsponge(ilci+2:nlci-1,:) )   
    281             spe2vr2(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) & 
    282                                            &        + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  )  
     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    ) )  
    283330         ENDIF 
    284331 
    285332         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    286             DO jj = 2, ispongearea 
    287                zlocalviscsponge(:,jj) = visc_dyn * ( ispongearea-jj ) * z1spongearea 
    288             ENDDO 
    289             spe1ur2(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) & 
    290                                              &      + zlocalviscsponge(2:jpi,2:ispongearea) )  
    291             spe2vr2(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     & 
    292                                              &      + zlocalviscsponge(:,3:ispongearea)     ) 
     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  ) ) 
    293339         ENDIF 
    294340 
    295341         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    296             DO jj = ilcj+1,nlcj-1 
    297                zlocalviscsponge(:,jj) = visc_dyn * (jj - (ilcj+1) ) * z1spongearea 
    298             ENDDO 
    299             spe1ur2(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) & 
    300                                            &         + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) )  
    301             spe2vr2(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) & 
    302                                            &         + zlocalviscsponge(:,ilcj+2:nlcj-1)     ) 
     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  ) ) 
    303348         ENDIF 
    304349         spongedoneU = .TRUE. 
    … …  
    306351      ENDIF 
    307352      ! 
    308       CALL wrk_dealloc( jpi, jpj, zlocalviscsponge ) 
     353      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp ) 
    309354      ! 
    310355#endif 
Note: See TracChangeset for help on using the changeset viewer.