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 12340 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_mus.F90 – NEMO

Ignore:
Timestamp:
2020-01-27T15:31:53+01:00 (4 years ago)
Author:
acc
Message:

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_mus.F90

    r12193 r12340  
    4747   !! * Substitutions 
    4848#  include "vectopt_loop_substitute.h90" 
     49#  include "do_loop_substitute.h90" 
    4950   !!---------------------------------------------------------------------- 
    5051   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    131132         zwx(:,:,jpk) = 0._wp                   ! bottom values 
    132133         zwy(:,:,jpk) = 0._wp   
    133          DO jk = 1, jpkm1                       ! interior values 
    134             DO jj = 1, jpjm1       
    135                DO ji = 1, fs_jpim1   ! vector opt. 
    136                   zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
    137                   zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
    138                END DO 
    139            END DO 
    140          END DO 
     134         DO_3D_10_10( 1, jpkm1 ) 
     135            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     136            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
     137         END_3D 
    141138         ! lateral boundary conditions   (changed sign) 
    142139         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. ) 
    … …  
    144141         zslpx(:,:,jpk) = 0._wp                 ! bottom values 
    145142         zslpy(:,:,jpk) = 0._wp 
    146          DO jk = 1, jpkm1                       ! interior values 
    147             DO jj = 2, jpj 
    148                DO ji = fs_2, jpi   ! vector opt. 
    149                   zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
    150                      &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
    151                   zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   & 
    152                      &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
    153                END DO 
    154             END DO 
    155          END DO 
    156          ! 
    157          DO jk = 1, jpkm1                 !-- Slopes limitation 
    158             DO jj = 2, jpj 
    159                DO ji = fs_2, jpi   ! vector opt. 
    160                   zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
    161                      &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
    162                      &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
    163                   zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
    164                      &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
    165                      &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
    166                END DO 
    167            END DO 
    168          END DO 
    169          ! 
    170          DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes 
    171             DO jj = 2, jpjm1 
    172                DO ji = fs_2, fs_jpim1   ! vector opt. 
    173                   ! MUSCL fluxes 
    174                   z0u = SIGN( 0.5, pU(ji,jj,jk) ) 
    175                   zalpha = 0.5 - z0u 
    176                   zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
    177                   zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
    178                   zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
    179                   zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    180                   ! 
    181                   z0v = SIGN( 0.5, pV(ji,jj,jk) ) 
    182                   zalpha = 0.5 - z0v 
    183                   zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    184                   zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
    185                   zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
    186                   zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    187                END DO 
    188             END DO 
    189          END DO 
     143         DO_3D_01_01( 1, jpkm1 ) 
     144            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
     145               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
     146            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   & 
     147               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
     148         END_3D 
     149         ! 
     150         DO_3D_01_01( 1, jpkm1 ) 
     151            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
     152               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
     153               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
     154            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
     155               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
     156               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
     157         END_3D 
     158         ! 
     159         DO_3D_00_00( 1, jpkm1 ) 
     160            ! MUSCL fluxes 
     161            z0u = SIGN( 0.5, pU(ji,jj,jk) ) 
     162            zalpha = 0.5 - z0u 
     163            zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     164            zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 
     165            zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk) 
     166            zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     167            ! 
     168            z0v = SIGN( 0.5, pV(ji,jj,jk) ) 
     169            zalpha = 0.5 - z0v 
     170            zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
     171            zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 
     172            zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk) 
     173            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     174         END_3D 
    190175         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign) 
    191176         ! 
    192          DO jk = 1, jpkm1                 !-- Tracer advective trend 
    193             DO jj = 2, jpjm1       
    194                DO ji = fs_2, fs_jpim1   ! vector opt. 
    195                   pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
    196                   &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
    197                   &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    198                END DO 
    199            END DO 
    200          END DO         
     177         DO_3D_00_00( 1, jpkm1 ) 
     178            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       & 
     179            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     & 
     180            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     181         END_3D 
    201182         !                                ! trend diagnostics 
    202183         IF( l_trd )  THEN 
    … …  
    219200         !                                !-- Slopes of tracer 
    220201         zslpx(:,:,1) = 0._wp                   ! surface values 
    221          DO jk = 2, jpkm1                       ! interior value 
    222             DO jj = 1, jpj 
    223                DO ji = 1, jpi 
    224                   zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
    225                      &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
    226                END DO 
    227             END DO 
    228          END DO 
    229          DO jk = 2, jpkm1                 !-- Slopes limitation 
    230             DO jj = 1, jpj                      ! interior values 
    231                DO ji = 1, jpi 
    232                   zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
    233                      &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
    234                      &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  ) 
    235                END DO 
    236             END DO 
    237          END DO 
    238          DO jk = 1, jpk-2                 !-- vertical advective flux 
    239             DO jj = 2, jpjm1       
    240                DO ji = fs_2, fs_jpim1   ! vector opt. 
    241                   z0w = SIGN( 0.5, pW(ji,jj,jk+1) ) 
    242                   zalpha = 0.5 + z0w 
    243                   zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm) 
    244                   zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
    245                   zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
    246                   zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
    247                END DO  
    248             END DO 
    249          END DO 
     202         DO_3D_11_11( 2, jpkm1 ) 
     203            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  & 
     204               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  ) 
     205         END_3D 
     206         DO_3D_11_11( 2, jpkm1 ) 
     207            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     208               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
     209               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  ) 
     210         END_3D 
     211         DO_3D_00_00( 1, jpk-2 ) 
     212            z0w = SIGN( 0.5, pW(ji,jj,jk+1) ) 
     213            zalpha = 0.5 + z0w 
     214            zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm) 
     215            zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 
     216            zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  ) 
     217            zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
     218         END_3D 
    250219         IF( ln_linssh ) THEN                   ! top values, linear free surface only 
    251220            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean) 
    252                DO jj = 1, jpj 
    253                   DO ji = 1, jpi 
    254                      zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) 
    255                   END DO 
    256                END DO    
     221               DO_2D_11_11 
     222                  zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) 
     223               END_2D 
    257224            ELSE                                      ! no cavities: only at the ocean surface 
    258225               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 
    … …  
    260227         ENDIF 
    261228         ! 
    262          DO jk = 1, jpkm1                 !-- vertical advective trend 
    263             DO jj = 2, jpjm1       
    264                DO ji = fs_2, fs_jpim1   ! vector opt. 
    265                   pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    266                END DO 
    267             END DO 
    268          END DO 
     229         DO_3D_00_00( 1, jpkm1 ) 
     230            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     231         END_3D 
    269232         !                                ! send trends for diagnostic 
    270233         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) ) 
Note: See TracChangeset for help on using the changeset viewer.