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 13295 for NEMO/trunk/src/OCE/TRA/traadv_fct.F90 – NEMO

Ignore:
Timestamp:
2020-07-10T20:24:21+02:00 (4 years ago)
Author:
acc
Message:

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/OCE/TRA/traadv_fct.F90

    r13286 r13295  
    139139      IF( ll_zAimp ) THEN 
    140140         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
    141          DO_3D_00_00( 1, jpkm1 ) 
     141         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    142142            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   & 
    143143            &                               / e3t(ji,jj,jk,Krhs) 
    … …  
    151151         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
    152152         !                    !* upstream tracer flux in the i and j direction  
    153          DO_3D_10_10( 1, jpkm1 ) 
     153         DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    154154            ! upstream scheme 
    155155            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) 
    … …  
    161161         END_3D 
    162162         !                    !* upstream tracer flux in the k direction *! 
    163          DO_3D_11_11( 2, jpkm1 ) 
     163         DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
    164164            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 
    165165            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 
    … …  
    168168         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked) 
    169169            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
    170                DO_2D_11_11 
     170               DO_2D( 1, 1, 1, 1 ) 
    171171                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface  
    172172               END_2D 
    173173            ELSE                             ! no cavities: only at the ocean surface 
    174                DO_2D_11_11 
     174               DO_2D( 1, 1, 1, 1 ) 
    175175                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 
    176176               END_2D 
    … …  
    178178         ENDIF 
    179179         !                
    180          DO_3D_00_00( 1, jpkm1 ) 
     180         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    181181            !                             ! total intermediate advective trends 
    182182            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    … …  
    194194            ! 
    195195            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
    196             DO_3D_00_00( 2, jpkm1 ) 
     196            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    197197               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    198198               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    … …  
    200200               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
    201201            END_3D 
    202             DO_3D_00_00( 1, jpkm1 ) 
     202            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    203203               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    204204                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    … …  
    218218         ! 
    219219         CASE(  2  )                   !- 2nd order centered 
    220             DO_3D_10_10( 1, jpkm1 ) 
     220            DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    221221               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk) 
    222222               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk) 
    … …  
    227227            zltv(:,:,jpk) = 0._wp 
    228228            DO jk = 1, jpkm1                 ! Laplacian 
    229                DO_2D_10_10 
     229               DO_2D( 1, 0, 1, 0 ) 
    230230                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
    231231                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
    232232               END_2D 
    233                DO_2D_00_00 
     233               DO_2D( 0, 0, 0, 0 ) 
    234234                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
    235235                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
    … …  
    238238            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
    239239            ! 
    240             DO_3D_10_10( 1, jpkm1 ) 
     240            DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    241241               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points 
    242242               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
    … …  
    249249            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
    250250            ztv(:,:,jpk) = 0._wp 
    251             DO_3D_10_10( 1, jpkm1 ) 
     251            DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    252252               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
    253253               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
    … …  
    255255            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
    256256            ! 
    257             DO_3D_00_00( 1, jpkm1 ) 
     257            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    258258               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2) 
    259259               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
    … …  
    271271         ! 
    272272         CASE(  2  )                   !- 2nd order centered 
    273             DO_3D_00_00( 2, jpkm1 ) 
     273            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    274274               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
    275275                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
    … …  
    278278         CASE(  4  )                   !- 4th order COMPACT 
    279279            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
    280             DO_3D_00_00( 2, jpkm1 ) 
     280            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    281281               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    282282            END_3D 
    … …  
    288288         !          
    289289         IF ( ll_zAimp ) THEN 
    290             DO_3D_00_00( 1, jpkm1 ) 
     290            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    291291               !                             ! total intermediate advective trends 
    292292               ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    … …  
    298298            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
    299299            ! 
    300             DO_3D_00_00( 2, jpkm1 ) 
     300            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    301301               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    302302               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    … …  
    313313         !        !==  final trend with corrected fluxes  ==! 
    314314         ! 
    315          DO_3D_00_00( 1, jpkm1 ) 
     315         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    316316            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    317317               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    … …  
    324324            ! 
    325325            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 
    326             DO_3D_00_00( 2, jpkm1 ) 
     326            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    327327               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    328328               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    … …  
    330330               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
    331331            END_3D 
    332             DO_3D_00_00( 1, jpkm1 ) 
     332            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    333333               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    334334                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    … …  
    409409      DO jk = 1, jpkm1 
    410410         ikm1 = MAX(jk-1,1) 
    411          DO_2D_00_00 
     411         DO_2D( 0, 0, 0, 0 ) 
    412412 
    413413            ! search maximum in neighbourhood 
    … …  
    443443      ! 3. monotonic flux in the i & j direction (paa & pbb) 
    444444      ! ---------------------------------------- 
    445       DO_3D_00_00( 1, jpkm1 ) 
     445      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    446446         zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    447447         zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
    … …  
    481481      !!---------------------------------------------------------------------- 
    482482       
    483       DO_3D_11_11( 3, jpkm1 ) 
     483      DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
    484484         zwd (ji,jj,jk) = 4._wp 
    485485         zwi (ji,jj,jk) = 1._wp 
    … …  
    496496      ! 
    497497      jk = 2                                          ! Switch to second order centered at top 
    498       DO_2D_11_11 
     498      DO_2D( 1, 1, 1, 1 ) 
    499499         zwd (ji,jj,jk) = 1._wp 
    500500         zwi (ji,jj,jk) = 0._wp 
    … …  
    504504      ! 
    505505      !                       !==  tridiagonal solve  ==! 
    506       DO_2D_11_11 
     506      DO_2D( 1, 1, 1, 1 ) 
    507507         zwt(ji,jj,2) = zwd(ji,jj,2) 
    508508      END_2D 
    509       DO_3D_11_11( 3, jpkm1 ) 
     509      DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
    510510         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    511511      END_3D 
    512512      ! 
    513       DO_2D_11_11 
     513      DO_2D( 1, 1, 1, 1 ) 
    514514         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    515515      END_2D 
    516       DO_3D_11_11( 3, jpkm1 ) 
     516      DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
    517517         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    518518      END_3D 
    519519 
    520       DO_2D_11_11 
     520      DO_2D( 1, 1, 1, 1 ) 
    521521         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    522522      END_2D 
    523       DO_3DS_11_11( jpk-2, 2, -1 ) 
     523      DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 
    524524         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    525525      END_3D 
    … …  
    546546      !                      !==  build the three diagonal matrix & the RHS  ==! 
    547547      ! 
    548       DO_3D_00_00( 3, jpkm1 ) 
     548      DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
    549549         zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal 
    550550         zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal 
    … …  
    565565      END IF 
    566566      ! 
    567       DO_2D_00_00 
     567      DO_2D( 0, 0, 0, 0 ) 
    568568         ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point 
    569569         ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point 
    … …  
    582582      !                       !==  tridiagonal solver  ==! 
    583583      ! 
    584       DO_2D_00_00 
     584      DO_2D( 0, 0, 0, 0 ) 
    585585         zwt(ji,jj,2) = zwd(ji,jj,2) 
    586586      END_2D 
    587       DO_3D_00_00( 3, jpkm1 ) 
     587      DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
    588588         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    589589      END_3D 
    590590      ! 
    591       DO_2D_00_00 
     591      DO_2D( 0, 0, 0, 0 ) 
    592592         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    593593      END_2D 
    594       DO_3D_00_00( 3, jpkm1 ) 
     594      DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 
    595595         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    596596      END_3D 
    597597 
    598       DO_2D_00_00 
     598      DO_2D( 0, 0, 0, 0 ) 
    599599         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    600600      END_2D 
    601       DO_3DS_00_00( jpk-2, 2, -1 ) 
     601      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 
    602602         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    603603      END_3D 
    … …  
    638638      kstart =  1  + klev 
    639639      ! 
    640       DO_2D_00_00 
     640      DO_2D( 0, 0, 0, 0 ) 
    641641         zwt(ji,jj,kstart) = pD(ji,jj,kstart) 
    642642      END_2D 
    643       DO_3D_00_00( kstart+1, jpkm1 ) 
     643      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 
    644644         zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    645645      END_3D 
    646646      ! 
    647       DO_2D_00_00 
     647      DO_2D( 0, 0, 0, 0 ) 
    648648         pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 
    649649      END_2D 
    650       DO_3D_00_00( kstart+1, jpkm1 ) 
     650      DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 
    651651         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    652652      END_3D 
    653653 
    654       DO_2D_00_00 
     654      DO_2D( 0, 0, 0, 0 ) 
    655655         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    656656      END_2D 
    657       DO_3DS_00_00( jpk-2, kstart, -1 ) 
     657      DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 ) 
    658658         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    659659      END_3D 
Note: See TracChangeset for help on using the changeset viewer.