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/SBC/sbcblk_phy.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/SBC/sbcblk_phy.F90

    r12182 r12340  
    100100   PUBLIC bulk_formula 
    101101 
     102   !! * Substitutions 
     103#  include "do_loop_substitute.h90" 
    102104   !!---------------------------------------------------------------------- 
    103105   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    179181      !!---------------------------------------------------------------------------------- 
    180182      ! 
    181       DO jj = 1, jpj 
    182          DO ji = 1, jpi 
    183             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    184             ztc2 = ztc*ztc 
    185             visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
    186          END DO 
    187       END DO 
     183      DO_2D_11_11 
     184         ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
     185         ztc2 = ztc*ztc 
     186         visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
     187      END_2D 
    188188      ! 
    189189   END FUNCTION visc_air 
    … …  
    270270      INTEGER  ::   ji, jj         ! dummy loop indices 
    271271      !!---------------------------------------------------------------------------------- 
    272       DO jj = 1, jpj 
    273          DO ji = 1, jpi 
    274             gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 
    275          END DO 
    276       END DO 
     272      DO_2D_11_11 
     273         gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 
     274      END_2D 
    277275   END FUNCTION gamma_moist_vctr 
    278276 
    … …  
    317315      !!------------------------------------------------------------------- 
    318316      ! 
    319       DO jj = 1, jpj 
    320          DO ji = 1, jpi 
    321             ! 
    322             zqa = (1._wp + rctv0*pqa(ji,jj)) 
    323             ! 
    324             ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
    325             !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
    326             !                      or 
    327             !  b/  -u* [ theta*              + 0.61 theta q* ] 
    328             ! 
    329             One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
    330                &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
    331             ! 
    332          END DO 
    333       END DO 
     317      DO_2D_11_11 
     318         ! 
     319         zqa = (1._wp + rctv0*pqa(ji,jj)) 
     320         ! 
     321         ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 
     322         !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 
     323         !                      or 
     324         !  b/  -u* [ theta*              + 0.61 theta q* ] 
     325         ! 
     326         One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 
     327            &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 
     328         ! 
     329      END_2D 
    334330      ! 
    335331      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 
    … …  
    355351      !!------------------------------------------------------------------- 
    356352      ! 
    357       DO jj = 1, jpj 
    358          DO ji = 1, jpi 
    359             ! 
    360             zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer... 
    361             zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
    362             zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
    363             zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer 
    364             ! 
    365             zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!) 
    366             ! 
    367             zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature" 
    368             ! 
    369             ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer 
    370             ! 
    371             Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk 
    372             ! 
    373          END DO 
    374       END DO 
     353      DO_2D_11_11 
     354         ! 
     355         zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer... 
     356         zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
     357         zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta,        zqa)*pz ) ! ~ mean absolute temperature of air within the layer 
     358         zgamma =  gamma_moist(zta, zqa)                                              ! Adiabatic lapse-rate for moist air within the layer 
     359         ! 
     360         zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!) 
     361         ! 
     362         zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature" 
     363         ! 
     364         ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) )  ! ~ mean absolute virtual temp. within the layer 
     365         ! 
     366         Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) )                            ! the usual definition of Ri_bulk 
     367         ! 
     368      END_2D 
    375369   END FUNCTION Ri_bulk 
    376370 
    … …  
    454448      !!---------------------------------------------------------------------------------- 
    455449      ! 
    456       DO jj = 1, jpj 
    457          DO ji = 1, jpi 
    458             ! 
    459             ze_sat =  e_sat_sclr( ptak(ji,jj) ) 
    460             ! 
    461             q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 
    462             ! 
    463          END DO 
    464       END DO 
     450      DO_2D_11_11 
     451         ! 
     452         ze_sat =  e_sat_sclr( ptak(ji,jj) ) 
     453         ! 
     454         q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 
     455         ! 
     456      END_2D 
    465457      ! 
    466458   END FUNCTION q_sat 
    … …  
    481473      !!---------------------------------------------------------------------------------- 
    482474      ! 
    483       DO jj = 1, jpj 
    484          DO ji = 1, jpi 
    485             ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
    486             q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
    487          END DO 
    488       END DO 
     475      DO_2D_11_11 
     476         ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
     477         q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
     478      END_2D 
    489479      ! 
    490480   END FUNCTION q_air_rh 
    … …  
    521511      INTEGER  ::   ji, jj     ! dummy loop indices 
    522512      !!---------------------------------------------------------------------------------- 
    523       DO jj = 1, jpj 
    524          DO ji = 1, jpi 
    525  
    526             zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
    527             zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
    528             zz0 = pust(ji,jj)/pUb(ji,jj) 
    529             zCd = zz0*zz0 
    530             zCh = zz0*ptst(ji,jj)/zdt 
    531             zCe = zz0*pqst(ji,jj)/zdq 
    532  
    533             CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
    534                &              zCd, zCh, zCe,                                        & 
    535                &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 & 
    536                &              pTau(ji,jj), zQsen, zQlat ) 
    537  
    538             zTs2  = pTs(ji,jj)*pTs(ji,jj) 
    539             zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
    540  
    541             pQns(ji,jj) = zQlat + zQsen + zQlw 
    542  
    543             IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
    544          END DO 
    545       END DO 
     513      DO_2D_11_11 
     514 
     515         zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
     516         zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
     517         zz0 = pust(ji,jj)/pUb(ji,jj) 
     518         zCd = zz0*zz0 
     519         zCh = zz0*ptst(ji,jj)/zdt 
     520         zCe = zz0*pqst(ji,jj)/zdq 
     521 
     522         CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
     523            &              zCd, zCh, zCe,                                        & 
     524            &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 & 
     525            &              pTau(ji,jj), zQsen, zQlat ) 
     526 
     527         zTs2  = pTs(ji,jj)*pTs(ji,jj) 
     528         zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
     529 
     530         pQns(ji,jj) = zQlat + zQsen + zQlw 
     531 
     532         IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
     533      END_2D 
    546534   END SUBROUTINE UPDATE_QNSOL_TAU 
    547535 
    … …  
    574562      INTEGER  :: ji, jj, jq     ! dummy loop indices 
    575563      !!---------------------------------------------------------------------------------- 
    576       DO jj = 1, jpj 
    577          DO ji = 1, jpi 
    578  
    579             !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
    580             ztaa = pTa(ji,jj) ! first guess... 
    581             DO jq = 1, 4 
    582                zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 
    583                ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder... 
    584             END DO 
    585             zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 
    586             zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
    587  
    588             zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10 
    589  
    590             pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 
    591  
    592             zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) 
    593             pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 
    594             pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 
    595  
    596             IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
    597             IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
    598  
     564      DO_2D_11_11 
     565 
     566         !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
     567         ztaa = pTa(ji,jj) ! first guess... 
     568         DO jq = 1, 4 
     569            zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 
     570            ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder... 
    599571         END DO 
    600       END DO 
     572         zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 
     573         zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
     574 
     575         zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10 
     576 
     577         pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 
     578 
     579         zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) 
     580         pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 
     581         pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 
     582 
     583         IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
     584         IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
     585 
     586      END_2D 
    601587   END SUBROUTINE BULK_FORMULA_VCTR 
    602588 
Note: See TracChangeset for help on using the changeset viewer.