Changeset 785


Ignore:
Timestamp:
2008-01-08T15:33:17+01:00 (13 years ago)
Author:
rblod
Message:

Improve traadv_tvd performance, see ticket #46

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r719 r785  
    7575      INTEGER  ::   ji, jj, jk              ! dummy loop indices 
    7676      REAL(wp) ::                        &  ! temporary scalar 
    77          ztai, ztaj, ztak,               &  !    "         "    
    78          zsai, zsaj, zsak,               &  !    "         "    
     77         ztat, zsat,                     &  !    "         "    
    7978         z_hdivn_x, z_hdivn_y, z_hdivn 
    8079      REAL(wp) ::   & 
     
    158157      ! total advective trend 
    159158      DO jk = 1, jpkm1 
     159         z2dtt = z2 * rdttra(jk) 
    160160         DO jj = 2, jpjm1 
    161161            DO ji = fs_2, fs_jpim1   ! vector opt. 
    162162               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    163                ! i- j- horizontal & k- vertical advective trends 
    164                ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  ) ) * zbtr 
    165                ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  ) ) * zbtr 
    166                ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr 
    167                zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  ) ) * zbtr 
    168                zsaj = - ( zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  ) ) * zbtr 
    169                zsak = - ( zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr 
    170163               ! total intermediate advective trends 
    171                zti(ji,jj,jk) = ztai + ztaj + ztak 
    172                zsi(ji,jj,jk) = zsai + zsaj + zsak 
     164               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   & 
     165                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   & 
     166                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr 
     167               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   &  
     168                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   & 
     169                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr 
     170               ! update and guess with monotonic sheme 
     171               ta(ji,jj,jk) =  ta(ji,jj,jk) + ztat 
     172               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsat 
     173               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ztat ) * tmask(ji,jj,jk) 
     174               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsat ) * tmask(ji,jj,jk) 
    173175            END DO 
    174176         END DO 
     
    222224         ! 
    223225      ENDIF 
    224  
    225       ! update and guess with monotonic sheme 
    226       DO jk = 1, jpkm1 
    227          z2dtt = z2 * rdttra(jk) 
    228          DO jj = 2, jpjm1 
    229             DO ji = fs_2, fs_jpim1   ! vector opt. 
    230                ta(ji,jj,jk) =  ta(ji,jj,jk) + zti(ji,jj,jk) 
    231                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsi(ji,jj,jk) 
    232                zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk) 
    233                zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsi(ji,jj,jk) ) * tmask(ji,jj,jk) 
    234             END DO 
    235          END DO 
    236       END DO 
    237226 
    238227      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
     
    290279            DO ji = fs_2, fs_jpim1   ! vector opt.   
    291280               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    292                ! i- j- horizontal & k- vertical advective trends 
    293                ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )) * zbtr 
    294                ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )) * zbtr 
    295                ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)) * zbtr 
    296                zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )) * zbtr 
    297                zsaj = - ( zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )) * zbtr 
    298                zsak = - ( zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1)) * zbtr 
    299  
     281               ! total advective trends 
     282               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   & 
     283                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   & 
     284                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr 
     285               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   & 
     286                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   & 
     287                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr 
    300288               ! add them to the general tracer trends 
    301                ta(ji,jj,jk) = ta(ji,jj,jk) + ztai + ztaj + ztak 
    302                sa(ji,jj,jk) = sa(ji,jj,jk) + zsai + zsaj + zsak 
     289               ta(ji,jj,jk) = ta(ji,jj,jk) + ztat 
     290               sa(ji,jj,jk) = sa(ji,jj,jk) + zsat 
    303291            END DO 
    304292         END DO 
     
    396384      !!       in-space based differencing for fluid 
    397385      !!---------------------------------------------------------------------- 
    398       REAL(wp), INTENT( in ) ::   prdt                               ! ??? 
     386      REAL(wp), INTENT( in ) ::   prdt   ! ??? 
    399387      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) ::   & 
    400388         pbef,                            & ! before field 
     
    407395      INTEGER ::   ikm1 
    408396      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo 
     397      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbup, zbdo 
    409398      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt 
     399      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv 
     400      REAL(wp) ::   zup, zdo 
    410401      !!---------------------------------------------------------------------- 
    411402 
    412403      zbig = 1.e+40 
    413404      zrtrn = 1.e-15 
    414       zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0   
     405      zbetup(:,:,jpk) = 0.e0   ;   zbetdo(:,:,jpk) = 0.e0 
     406 
    415407 
    416408      ! Search local extrema 
    417409      ! -------------------- 
    418       ! large negative value (-zbig) inside land 
    419       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    420       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 
    421       ! search maximum in neighbourhood 
     410      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
     411      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   & 
     412         &        paft * tmask - zbig * ( 1.e0 - tmask )  ) 
     413      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   & 
     414         &        paft * tmask + zbig * ( 1.e0 - tmask )  ) 
     415 
    422416      DO jk = 1, jpkm1 
    423417         ikm1 = MAX(jk-1,1) 
    424          DO jj = 2, jpjm1 
    425             DO ji = fs_2, fs_jpim1   ! vector opt. 
    426                zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
    427                   &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   & 
    428                   &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   & 
    429                   &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   & 
    430                   &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   & 
    431                   &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
    432                   &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
    433             END DO 
    434          END DO 
    435       END DO 
    436       ! large positive value (+zbig) inside land 
    437       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    438       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 
    439       ! search minimum in neighbourhood 
    440       DO jk = 1, jpkm1 
    441          ikm1 = MAX(jk-1,1) 
    442          DO jj = 2, jpjm1 
    443             DO ji = fs_2, fs_jpim1   ! vector opt. 
    444                zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   & 
    445                   &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   & 
    446                   &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   & 
    447                   &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   & 
    448                   &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   & 
    449                   &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   & 
    450                   &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  ) 
    451             END DO 
    452          END DO 
    453       END DO 
    454  
    455       ! restore masked values to zero 
    456       pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 
    457       paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 
    458   
    459  
    460       ! 2. Positive and negative part of fluxes and beta terms 
    461       ! ------------------------------------------------------ 
    462  
    463       DO jk = 1, jpkm1 
    464418         z2dtt = prdt * rdttra(jk) 
    465419         DO jj = 2, jpjm1 
    466420            DO ji = fs_2, fs_jpim1   ! vector opt. 
    467                ! positive & negative part of the flux 
     421 
     422               ! search maximum in neighbourhood 
     423               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
     424                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
     425                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   & 
     426                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  ) 
     427 
     428               ! search minimum in neighbourhood 
     429               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
     430                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
     431                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   & 
     432                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  ) 
     433 
     434               ! positive part of the flux 
    468435               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
    469436                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   & 
    470437                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
     438 
     439               ! negative part of the flux 
    471440               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
    472441                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   & 
    473442                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
     443 
    474444               ! up & down beta terms 
    475445               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 
    476                zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 
    477                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt 
     446               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
     447               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
    478448            END DO 
    479449         END DO 
     
    490460         DO jj = 2, jpjm1 
    491461            DO ji = fs_2, fs_jpim1   ! vector opt. 
    492                za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    493                zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
    494                zc = 0.5 * ( 1.e0 + SIGN( 1.e0, paa(ji,jj,jk) ) ) 
    495                paa(ji,jj,jk) = paa(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) 
    496  
    497                za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
    498                zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
    499                zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pbb(ji,jj,jk) ) ) 
    500                pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) 
    501             END DO 
    502          END DO 
    503       END DO 
    504  
     462               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     463               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     464               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
     465               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 
     466 
     467               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     468               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     469               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
     470               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 
    505471 
    506472      ! monotonic flux in the k direction, i.e. pcc 
    507473      ! ------------------------------------------- 
    508       DO jk = 2, jpkm1 
    509          DO jj = 2, jpjm1 
    510             DO ji = fs_2, fs_jpim1   ! vector opt. 
    511  
    512                za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) ) 
    513                zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) ) 
    514                zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) ) 
    515                pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) 
     474               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 
     475               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
     476               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 
     477               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb ) 
    516478            END DO 
    517479         END DO 
     
    521483      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign 
    522484      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign 
    523       CALL lbc_lnk( pcc, 'W',  1. )      ! NO changed sign 
    524485      ! 
    525486   END SUBROUTINE nonosc 
Note: See TracChangeset for help on using the changeset viewer.