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/DIU/diu_bulk.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/DIU/diu_bulk.F90

    r11960 r12340  
    3636 
    3737   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step 
     38   !! * Substitutions 
     39#  include "do_loop_substitute.h90" 
    3840    
    3941   !!---------------------------------------------------------------------- 
    … …  
    128130      ! If not done already, calculate the solar fraction 
    129131      IF ( kt==nit000 ) THEN 
    130          DO jj = 1,jpj 
    131             DO ji = 1, jpi 
    132                IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) & 
    133                   &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) ) 
    134             END DO 
    135          END DO    
     132         DO_2D_11_11 
     133            IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) & 
     134               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) ) 
     135         END_2D 
    136136      ENDIF    
    137137 
    … …  
    199199      INTEGER :: ji,jj 
    200200                      
    201       DO jj = 1, jpj 
    202          DO ji = 1, jpi    
     201      DO_2D_11_11 
     202          
     203         ! Only calculate outside tmask 
     204         IF ( tmask(ji,jj,1) /= 1._wp ) THEN 
     205            t_imp(ji,jj) = 0._wp 
     206            CYCLE    
     207         END IF 
     208          
     209         IF (p_fvel(ji,jj) < pp_min_fvel) THEN 
     210            z_fvel = pp_min_fvel 
     211            WRITE(warn_string,*) "diurnal_sst_takaya step: "& 
     212            &//"friction velocity < minimum\n" & 
     213            &//"Setting friction velocity =",pp_min_fvel 
     214            CALL ctl_warn(warn_string) 
    203215             
    204             ! Only calculate outside tmask 
    205             IF ( tmask(ji,jj,1) /= 1._wp ) THEN 
    206                t_imp(ji,jj) = 0._wp 
    207                CYCLE    
    208             END IF 
    209              
    210             IF (p_fvel(ji,jj) < pp_min_fvel) THEN 
    211                z_fvel = pp_min_fvel 
    212                WRITE(warn_string,*) "diurnal_sst_takaya step: "& 
    213                &//"friction velocity < minimum\n" & 
    214                &//"Setting friction velocity =",pp_min_fvel 
    215                CALL ctl_warn(warn_string) 
    216                 
    217             ELSE 
    218                z_fvel = p_fvel(ji,jj) 
    219             ENDIF 
    220                   
    221             ! Calculate the Obukhov length 
    222             IF ( (z_fvel < pp_veltol ) .AND. & 
    223             &    (p_dsst(ji,jj) > 0._wp) ) THEN 
    224                z_olength =  z_fvel  / & 
    225                &     SQRT( p_dsst(ji,jj) * vkarmn * grav * & 
    226                &             pp_alpha / ( 5._wp * pthick(ji,jj) ) ) 
    227             ELSE 
    228                z_olength = & 
    229                &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / & 
    230                &   ( vkarmn * grav * pp_alpha *& 
    231                &     p_abflux(ji,jj) ) 
    232             ENDIF 
    233               
    234             ! Calculate the stability function           
    235             z_sigma = pthick(ji,jj) / z_olength 
    236             z_sigma2 = z_sigma * z_sigma 
    237    
    238             IF ( z_sigma >= 0. ) THEN 
    239                z_stabfunc = 1._wp + & 
    240                &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / & 
    241                &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * & 
    242                &    z_sigma2 ) ) 
    243             ELSE 
    244                z_stabfunc = 1._wp / & 
    245                &                   SQRT( 1._wp - 16._wp * z_sigma ) 
    246             ENDIF 
    247  
    248             ! Calculate the T increment 
    249             z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / & 
    250             & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) ) 
    251              
    252               
    253             z_term2 = -( ( pmu(ji,jj) + 1._wp) * & 
    254             &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / & 
    255             &      ( pthick(ji,jj) * z_stabfunc ) )      
     216         ELSE 
     217            z_fvel = p_fvel(ji,jj) 
     218         ENDIF 
     219               
     220         ! Calculate the Obukhov length 
     221         IF ( (z_fvel < pp_veltol ) .AND. & 
     222         &    (p_dsst(ji,jj) > 0._wp) ) THEN 
     223            z_olength =  z_fvel  / & 
     224            &     SQRT( p_dsst(ji,jj) * vkarmn * grav * & 
     225            &             pp_alpha / ( 5._wp * pthick(ji,jj) ) ) 
     226         ELSE 
     227            z_olength = & 
     228            &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / & 
     229            &   ( vkarmn * grav * pp_alpha *& 
     230            &     p_abflux(ji,jj) ) 
     231         ENDIF 
    256232           
    257             t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / & 
    258                            ( 1._wp - p_rdt * z_term2 ) 
    259  
    260          END DO 
    261       END DO 
     233         ! Calculate the stability function           
     234         z_sigma = pthick(ji,jj) / z_olength 
     235         z_sigma2 = z_sigma * z_sigma 
     236         IF ( z_sigma >= 0. ) THEN 
     237            z_stabfunc = 1._wp + & 
     238            &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / & 
     239            &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * & 
     240            &    z_sigma2 ) ) 
     241         ELSE 
     242            z_stabfunc = 1._wp / & 
     243            &                   SQRT( 1._wp - 16._wp * z_sigma ) 
     244         ENDIF 
     245 
     246         ! Calculate the T increment 
     247         z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / & 
     248         & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) ) 
     249          
     250           
     251         z_term2 = -( ( pmu(ji,jj) + 1._wp) * & 
     252         &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / & 
     253         &      ( pthick(ji,jj) * z_stabfunc ) )      
     254        
     255         t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / & 
     256                        ( 1._wp - p_rdt * z_term2 ) 
     257 
     258      END_2D 
    262259       
    263260      END FUNCTION t_imp 
Note: See TracChangeset for help on using the changeset viewer.