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/DIA/diahth.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/DIA/diahth.F90

    r12193 r12340  
    4040 
    4141 
     42   !! * Substitutions 
     43#  include "do_loop_substitute.h90" 
    4244   !!---------------------------------------------------------------------- 
    4345   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    127129            zdepinv(:,:) = 0._wp   
    128130            zmaxdzT(:,:) = 0._wp   
    129             DO jj = 1, jpj 
    130                DO ji = 1, jpi 
     131            DO_2D_11_11 
     132               zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
     133               hth     (ji,jj) = zztmp 
     134               zabs2   (ji,jj) = zztmp 
     135               ztm2    (ji,jj) = zztmp 
     136               zrho10_3(ji,jj) = zztmp 
     137               zpycn   (ji,jj) = zztmp 
     138            END_2D 
     139            IF( nla10 > 1 ) THEN  
     140               DO_2D_11_11 
    131141                  zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
    132                   hth     (ji,jj) = zztmp 
    133                   zabs2   (ji,jj) = zztmp 
    134                   ztm2    (ji,jj) = zztmp 
    135                   zrho10_3(ji,jj) = zztmp 
    136                   zpycn   (ji,jj) = zztmp 
    137                  END DO 
    138             END DO 
    139             IF( nla10 > 1 ) THEN  
    140                DO jj = 1, jpj 
    141                   DO ji = 1, jpi 
    142                      zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
    143                      zrho0_3(ji,jj) = zztmp 
    144                      zrho0_1(ji,jj) = zztmp 
    145                   END DO 
    146                END DO 
     142                  zrho0_3(ji,jj) = zztmp 
     143                  zrho0_1(ji,jj) = zztmp 
     144               END_2D 
    147145            ENDIF 
    148146       
    149147            ! Preliminary computation 
    150148            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    151             DO jj = 1, jpj 
    152                DO ji = 1, jpi 
    153                   IF( tmask(ji,jj,nla10) == 1. ) THEN 
    154                      zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
    155                         &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   & 
    156                         &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 
    157                      zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
    158                         &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 
    159                      zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm) 
    160                      zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 
    161                      zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    162                      zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    163                   ELSE 
    164                      zdelr(ji,jj) = 0._wp 
    165                   ENDIF 
    166                END DO 
    167             END DO 
     149            DO_2D_11_11 
     150               IF( tmask(ji,jj,nla10) == 1. ) THEN 
     151                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
     152                     &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   & 
     153                     &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 
     154                  zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
     155                     &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 
     156                  zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm) 
     157                  zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 
     158                  zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     159                  zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     160               ELSE 
     161                  zdelr(ji,jj) = 0._wp 
     162               ENDIF 
     163            END_2D 
    168164 
    169165            ! ------------------------------------------------------------- ! 
    … …  
    173169            ! MLD: rho = rho(1) + zrho1                                     ! 
    174170            ! ------------------------------------------------------------- ! 
    175             DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    176                DO jj = 1, jpj 
    177                   DO ji = 1, jpi 
    178                      ! 
    179                      zzdep = gdepw(ji,jj,jk,Kmm) 
    180                      zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) & 
    181                             & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    182                      zzdep = zzdep * tmask(ji,jj,1) 
    183  
    184                      IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
    185                          zmaxdzT(ji,jj) = zztmp    
    186                          hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    187                      ENDIF 
    188                 
    189                      IF( nla10 > 1 ) THEN  
    190                         zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1) 
    191                         IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03 
    192                         IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01 
    193                      ENDIF 
    194                   END DO 
    195                END DO 
    196             END DO 
     171            DO_3DS_11_11( jpkm1, 2, -1 ) 
     172               ! 
     173               zzdep = gdepw(ji,jj,jk,Kmm) 
     174               zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) & 
     175                      & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     176               zzdep = zzdep * tmask(ji,jj,1) 
     177 
     178               IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
     179                   zmaxdzT(ji,jj) = zztmp    
     180                   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     181               ENDIF 
     182          
     183               IF( nla10 > 1 ) THEN  
     184                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1) 
     185                  IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03 
     186                  IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01 
     187               ENDIF 
     188            END_3D 
    197189          
    198190            CALL iom_put( 'mlddzt', hth )            ! depth of the thermocline 
    … …  
    214206            ! depth of temperature inversion                                ! 
    215207            ! ------------------------------------------------------------- ! 
    216             DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    217                DO jj = 1, jpj 
    218                   DO ji = 1, jpi 
    219                      ! 
    220                      zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 
    221                      ! 
    222                      zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm)  ! - delta T(10m) 
    223                      IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    224                      IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    225                      zztmp = -zztmp                                          ! delta T(10m) 
    226                      IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    227                         ztinv(ji,jj) = zztmp    
    228                         zdepinv (ji,jj) = zzdep   ! max value and depth 
    229                      ENDIF 
    230  
    231                      zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    232                      IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    233                      IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    234                      ! 
    235                   END DO 
    236                END DO 
    237             END DO 
     208            DO_3DS_11_11( jpkm1, nlb10, -1 ) 
     209               ! 
     210               zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 
     211               ! 
     212               zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm)  ! - delta T(10m) 
     213               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     214               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     215               zztmp = -zztmp                                          ! delta T(10m) 
     216               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     217                  ztinv(ji,jj) = zztmp    
     218                  zdepinv (ji,jj) = zzdep   ! max value and depth 
     219               ENDIF 
     220 
     221               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     222               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     223               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     224               ! 
     225            END_3D 
    238226 
    239227            CALL iom_put( 'mld_dt02', zabs2    )   ! MLD abs(delta t) - 0.2 
    … …  
    316304      ! --------------------------------------- ! 
    317305      iktem(:,:) = 1 
    318       DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    319          DO jj = 1, jpj 
    320             DO ji = 1, jpi 
    321                zztmp = ts(ji,jj,jk,jp_tem,Kmm) 
    322                IF( zztmp >= ptem )   iktem(ji,jj) = jk 
    323             END DO 
    324          END DO 
    325       END DO 
     306      DO_3D_11_11( 1, jpkm1 ) 
     307         zztmp = ts(ji,jj,jk,jp_tem,Kmm) 
     308         IF( zztmp >= ptem )   iktem(ji,jj) = jk 
     309      END_3D 
    326310 
    327311      ! ------------------------------- ! 
    328312      !  Depth of ptem isotherm         ! 
    329313      ! ------------------------------- ! 
    330       DO jj = 1, jpj 
    331          DO ji = 1, jpi 
    332             ! 
    333             zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom 
    334             ! 
    335             iid = iktem(ji,jj) 
    336             IF( iid /= 1 ) THEN  
    337                 zztmp =     gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation 
    338                   &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   & 
    339                   &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   & 
    340                   &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) ) 
    341                pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    342             ELSE  
    343                pdept(ji,jj) = 0._wp 
    344             ENDIF 
    345          END DO 
    346       END DO 
     314      DO_2D_11_11 
     315         ! 
     316         zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom 
     317         ! 
     318         iid = iktem(ji,jj) 
     319         IF( iid /= 1 ) THEN  
     320             zztmp =     gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation 
     321               &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   & 
     322               &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   & 
     323               &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) ) 
     324            pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
     325         ELSE  
     326            pdept(ji,jj) = 0._wp 
     327         ENDIF 
     328      END_2D 
    347329      ! 
    348330   END SUBROUTINE dia_hth_dep 
    … …  
    368350      ! 
    369351      ilevel(:,:) = 1 
    370       DO jk = 2, jpkm1 
    371          DO jj = 1, jpj 
    372             DO ji = 1, jpi 
    373                IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 
    374                    ilevel(ji,jj) = jk 
    375                    zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 
    376                    phtc  (ji,jj) = phtc  (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk) 
    377                ENDIF 
    378             ENDDO 
    379          ENDDO 
    380       ENDDO 
    381       ! 
    382       DO jj = 1, jpj 
    383          DO ji = 1, jpi 
    384             ik = ilevel(ji,jj) 
    385             zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep 
    386             phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
    387                                                           * tmask(ji,jj,ik+1) 
    388          END DO 
    389       ENDDO 
     352      DO_3D_11_11( 2, jpkm1 ) 
     353         IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 
     354             ilevel(ji,jj) = jk 
     355             zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 
     356             phtc  (ji,jj) = phtc  (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk) 
     357         ENDIF 
     358      END_3D 
     359      ! 
     360      DO_2D_11_11 
     361         ik = ilevel(ji,jj) 
     362         zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep 
     363         phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
     364                                                       * tmask(ji,jj,ik+1) 
     365      END_2D 
    390366      ! 
    391367      ! 
Note: See TracChangeset for help on using the changeset viewer.