Changeset 2508 for LMDZ5/trunk


Ignore:
Timestamp:
May 7, 2016, 7:12:33 PM (9 years ago)
Author:
jyg
Message:

Optimization of cv3_yield (in cv3_routines.F90):
suppression of the triple vertical loops. The
optimized code is activated by:

ok_optim_yield=y

which changes numerical results.

Location:
LMDZ5/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cv3_routines.F90

    r2490 r2508  
    105105     ok_intermittent=.False.
    106106     CALL getin_p('ok_intermittent',ok_intermittent)
     107     ok_optim_yield=.False.
     108     CALL getin_p('ok_optim_yield',ok_optim_yield)
    107109     coef_peel=0.25
    108110     CALL getin_p('coef_peel',coef_peel)
    … …  
    126128    WRITE (*, *) 'tau_stop=', tau_stop
    127129    WRITE (*, *) 'ok_intermittent=', ok_intermittent
     130    WRITE (*, *) 'ok_optim_yield =', ok_optim_yield
    128131    WRITE (*, *) 'coef_peel=', coef_peel
    129132
    … …  
    29772980!
    29782981!local variables:
    2979       INTEGER i, k, il, n, j, num1
    2980       REAL rat, delti
    2981       REAL ax, bx, cx, dx, ex
    2982       REAL cpinv, rdcp, dpinv
    2983       REAL awat(nloc)
    2984       REAL lvcp(nloc, na), lfcp(nloc, na)                  ! , mke(nloc, na) ! unused . jyg
    2985       REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
     2982      INTEGER                                            :: i, k, il, n, j, num1
     2983      REAL                                               :: rat, delti
     2984      REAL                                               :: ax, bx, cx, dx, ex
     2985      REAL                                               :: cpinv, rdcp, dpinv
     2986      REAL, DIMENSION (nloc)                             ::  awat
     2987      REAL, DIMENSION (nloc, nd)                         :: lvcp, lfcp              ! , mke ! unused . jyg
     2988      REAL, DIMENSION (nloc)                             :: am, work, ad, amp1
    29862989!!      real up1(nloc), dn1(nloc)
    2987       REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
    2988       REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
    2989       REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
    2990       REAL th_wake(nloc, nd)
    2991       REAL alpha_qpos(nloc), alpha_qpos1(nloc)
    2992       REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
    2993       REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
    2994       REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
    2995       REAL qnk(nloc)
     2990      REAL, DIMENSION (nloc, nd, nd)                     :: up1, dn1
     2991!jyg<
     2992      REAL, DIMENSION (nloc, nd)                         :: up_to, up_from
     2993      REAL, DIMENSION (nloc, nd)                         :: dn_to, dn_from
     2994!>jyg
     2995      REAL, DIMENSION (nloc)                             :: asum, bsum, csum, dsum
     2996      REAL, DIMENSION (nloc)                             :: esum, fsum, gsum, hsum
     2997      REAL, DIMENSION (nloc, nd)                         :: th_wake
     2998      REAL, DIMENSION (nloc)                             :: alpha_qpos, alpha_qpos1
     2999      REAL, DIMENSION (nloc, nd)                         :: qcond, nqcond, wa           ! cld
     3000      REAL, DIMENSION (nloc, nd)                         :: siga, sax, mac              ! cld
     3001      REAL, DIMENSION (nloc)                             :: sument
     3002      REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
     3003      REAL, DIMENSION (nloc)                             :: qnk
    29963004      REAL sumdq !jyg
    29973005!
    … …  
    32333241! print*,'cv3_yield apres ft'
    32343242
     3243!jyg<
     3244!-----------------------------------------------------------
     3245           IF (ok_optim_yield) THEN                       !|
     3246!-----------------------------------------------------------
     3247!
     3248!***                                                      ***
     3249!***    Compute convective mass fluxes upwd and dnwd      ***
     3250
     3251upwd(:,:) = 0.
     3252up_to(:,:) = 0.
     3253up_from(:,:) = 0.
     3254dnwd(:,:) = 0.
     3255dn_to(:,:) = 0.
     3256dn_from(:,:) = 0.
     3257!
     3258! =================================================
     3259!              upward fluxes                      |
     3260! ------------------------------------------------
     3261DO i = 2, nl
     3262  DO il = 1, ncum
     3263    IF (i<=inb(il)) THEN
     3264      up_to(il,i) = m(il,i)
     3265    ENDIF
     3266  ENDDO
     3267  DO j = 1, i-1
     3268    DO il = 1, ncum
     3269      IF (i<=inb(il)) THEN
     3270        up_to(il,i) = up_to(il,i) + ment(il,j,i)
     3271      ENDIF
     3272    ENDDO
     3273  ENDDO
     3274ENDDO
     3275!
     3276DO i = 1, nl
     3277  DO il = 1, ncum
     3278    IF (i<=inb(il)) THEN
     3279      up_from(il,i) = cbmf(il)*wghti(il,i)
     3280    ENDIF
     3281  ENDDO
     3282ENDDO
     3283!!DO i = 2, nl
     3284!!  DO j = i+1, nl          !! Permuter les boucles i et j
     3285DO j = 3, nl
     3286  DO i = 2, j-1
     3287    DO il = 1, ncum
     3288      IF (j<=inb(il)) THEN
     3289        up_from(il,i) = up_from(il,i) + ment(il,i,j)
     3290      ENDIF
     3291    ENDDO
     3292  ENDDO
     3293ENDDO
     3294!
     3295! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer
     3296!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting
     3297!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)):
     3298!
     3299DO i = 2, nlp
     3300  DO il = 1, ncum
     3301    upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
     3302  ENDDO
     3303ENDDO
     3304!
     3305! =================================================
     3306!              downward fluxes                    |
     3307! ------------------------------------------------
     3308DO i = 1, nl
     3309  DO j = i+1, nl
     3310    DO il = 1, ncum
     3311      IF (j<=inb(il)) THEN
     3312        dn_to(il,i) = dn_to(il,i) + ment(il,j,i)
     3313      ENDIF
     3314    ENDDO
     3315  ENDDO
     3316ENDDO
     3317!
     3318!!DO i = 2, nl
     3319!!  DO j = 1, i-1          !! Permuter les boucles i et j
     3320DO j = 1, nl
     3321  DO i = j+1, nl
     3322    DO il = 1, ncum
     3323      IF (i<=inb(il)) THEN
     3324        dn_from(il,i) = dn_from(il,i) + ment(il,i,j)
     3325      ENDIF
     3326    ENDDO
     3327  ENDDO
     3328ENDDO
     3329!
     3330! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer
     3331!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts
     3332!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)):
     3333!
     3334DO i = nl-1, 1, -1
     3335  DO il = 1, ncum
     3336    dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
     3337  ENDDO
     3338ENDDO
     3339! =================================================
     3340!
     3341!-----------------------------------------------------------
     3342        ENDIF !(ok_optim_yield)                           !|
     3343!-----------------------------------------------------------
     3344!>jyg
     3345
    32353346! ***  calculate tendencies of potential temperature and mixing ratio  ***
    32363347! ***               at levels above the lowest level                   ***
    … …  
    32403351
    32413352
    3242   DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
     3353!jyg<
     3354!!  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
     3355  DO i = 2, nl
     3356!>jyg
    32433357
    32443358    num1 = 0
    … …  
    32483362    IF (num1<=0) GO TO 500
    32493363
     3364!
    32503365!jyg<
    3251 !!    CALL zilch(amp1, ncum)
    3252 !!    CALL zilch(ad, ncum)
     3366!-----------------------------------------------------------
     3367           IF (ok_optim_yield) THEN                       !|
     3368!-----------------------------------------------------------
     3369DO il = 1, ncum
     3370   amp1(il) = upwd(il,i+1)
     3371   ad(il) = dnwd(il,i)
     3372ENDDO
     3373!-----------------------------------------------------------
     3374        ELSE !(ok_optim_yield)                            !|
     3375!-----------------------------------------------------------
     3376!>jyg
    32533377    DO il = 1,ncum
    32543378      amp1(il) = 0.
    32553379      ad(il) = 0.
    32563380    ENDDO
    3257 !>jyg
    32583381
    32593382    DO k = 1, nl + 1
    … …  
    32723395    END DO
    32733396
    3274     DO k = 1, i
    3275       DO j = i + 1, nl + 1
     3397    DO j = i + 1, nl + 1         
     3398       DO k = 1, i
     3399          !yor! reverted j and k loops
     3400          DO il = 1, ncum
     3401!yor!        IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
     3402             IF (j<=(inb(il)+1)) THEN 
     3403                amp1(il) = amp1(il) + ment(il, k, j)
     3404             END IF
     3405          END DO
     3406       END DO
     3407    END DO
     3408
     3409    DO k = 1, i - 1
     3410!jyg<
     3411!!      DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
     3412      DO j = i, nl
     3413!>jyg
    32763414        DO il = 1, ncum
    3277           IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN
    3278             amp1(il) = amp1(il) + ment(il, k, j)
    3279           END IF
    3280         END DO
    3281       END DO
    3282     END DO
    3283 
    3284     DO k = 1, i - 1
    3285       DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
    3286         DO il = 1, ncum
    3287           IF (i<=inb(il) .AND. j<=inb(il)) THEN
     3415!yor!        IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
     3416             IF (j<=inb(il)) THEN   
    32883417            ad(il) = ad(il) + ment(il, j, k)
    32893418          END IF
    … …  
    32913420      END DO
    32923421    END DO
     3422!
     3423!-----------------------------------------------------------
     3424        ENDIF !(ok_optim_yield)                           !|
     3425!-----------------------------------------------------------
     3426!
     3427!!   print *,'yield, i, amp1, ad', i, amp1(1), ad(1)
    32933428
    32943429    DO il = 1, ncum
    … …  
    34383573!AC!      enddo
    34393574
    3440     DO k = i, nl + 1
     3575!jyg<
     3576!!    DO k = i, nl + 1
     3577    DO k = i, nl
     3578!>jyg
    34413579
    34423580      IF (iflag_mix/=0) THEN
    … …  
    37413879    END DO
    37423880  END DO
     3881!jyg<
     3882!-----------------------------------------------------------
     3883           IF (ok_optim_yield) THEN                       !|
     3884!-----------------------------------------------------------
    37433885  DO i = 1, nl
    3744     DO j = 1, nl
     3886    DO il = 1, ncum
     3887      IF (iflag(il)<=1) THEN
     3888        upwd(il, i) = upwd(il, i)/alpha_qpos(il)
     3889        dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
     3890      END IF
     3891    END DO
     3892  END DO
     3893!-----------------------------------------------------------
     3894        ENDIF !(ok_optim_yield)                           !|
     3895!-----------------------------------------------------------
     3896!>jyg
     3897  DO j = 1, nl !yor! inverted i and j loops
     3898     DO i = 1, nl
    37453899      DO il = 1, ncum
    37463900        IF (iflag(il)<=1) THEN
    … …  
    37763930  DO i = 1, nl
    37773931    DO il = 1, ncum
    3778       upwd(il, i) = 0.0
    3779       dnwd(il, i) = 0.0
    3780     END DO
    3781   END DO
    3782 
    3783   DO i = 1, nl
    3784     DO il = 1, ncum
    37853932      dnwd0(il, i) = -mp(il, i)
    37863933    END DO
    … …  
    37953942
    37963943
     3944!jyg<
     3945!-----------------------------------------------------------
     3946           IF (.NOT.ok_optim_yield) THEN                  !|
     3947!-----------------------------------------------------------
    37973948  DO i = 1, nl
    37983949    DO il = 1, ncum
    3799       IF (i>=icb(il) .AND. i<=inb(il)) THEN
    3800         upwd(il, i) = 0.0
    3801         dnwd(il, i) = 0.0
    3802       END IF
    3803     END DO
    3804   END DO
     3950      upwd(il, i) = 0.0
     3951      dnwd(il, i) = 0.0
     3952    END DO
     3953  END DO
     3954
     3955!!  DO i = 1, nl                                           ! useless; jyg
     3956!!    DO il = 1, ncum                                      ! useless; jyg
     3957!!      IF (i>=icb(il) .AND. i<=inb(il)) THEN              ! useless; jyg
     3958!!        upwd(il, i) = 0.0                                ! useless; jyg
     3959!!        dnwd(il, i) = 0.0                                ! useless; jyg
     3960!!      END IF                                             ! useless; jyg
     3961!!    END DO                                               ! useless; jyg
     3962!!  END DO                                                 ! useless; jyg
    38053963
    38063964  DO i = 1, nl
    … …  
    38133971  END DO
    38143972
     3973!yor! commented original
     3974!  DO i = 1, nl
     3975!    DO k = i, nl
     3976!      DO n = 1, i - 1
     3977!        DO il = 1, ncum
     3978!          IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
     3979!            up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
     3980!            dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
     3981!          END IF
     3982!        END DO
     3983!      END DO
     3984!    END DO
     3985!  END DO
     3986!yor! replaced with
    38153987  DO i = 1, nl
    38163988    DO k = i, nl
    38173989      DO n = 1, i - 1
    38183990        DO il = 1, ncum
    3819           IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
    3820             up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
    3821             dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
     3991          IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
     3992             up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
    38223993          END IF
    38233994        END DO
    … …  
    38253996    END DO
    38263997  END DO
     3998  DO i = 1, nl
     3999    DO n = 1, i - 1
     4000      DO k = i, nl
     4001        DO il = 1, ncum
     4002          IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor !  i always <= k
     4003             dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
     4004          END IF
     4005        END DO
     4006      END DO
     4007    END DO
     4008  END DO
     4009!yor! end replace
    38274010
    38284011  DO i = 1, nl
    … …  
    38754058!!!!
    38764059!!!!      ENDDO
     4060!-----------------------------------------------------------
     4061        ENDIF !(.NOT.ok_optim_yield)                      !|
     4062!-----------------------------------------------------------
     4063!>jyg
    38774064
    38784065! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  • LMDZ5/trunk/libf/phylmd/cv3param.h

    r2458 r2508  
    77!------------------------------------------------------------
    88
     9      logical ok_optim_yield
    910      logical ok_convstop
    1011      logical ok_intermittent
    … …  
    3637                      ,flag_wb &
    3738                      ,noff, minorig, nl, nlp, nlm  &
    38                       ,ok_convstop, ok_intermittent
     39                      ,ok_convstop, ok_intermittent &
     40                      ,ok_optim_yield
    3941!$OMP THREADPRIVATE(/cv3param/)
    4042
Note: See TracChangeset for help on using the changeset viewer.