Ignore:
Timestamp:
Jun 11, 2020, 11:09:41 AM (4 years ago)
Author:
adurocher
Message:

cv3p_mixing_new : merge some loops

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90

    r3710 r3711  
    131131
    132132    INTEGER, SAVE                                       :: igout = 1
    133     !$omp THREADPRIVATE(igout)
    134 
    135     ! --   Mixing probability distribution functions
     133!$omp THREADPRIVATE(igout)
     134
     135! --   Mixing probability distribution functions
    136136
    137137    REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
     
    149149
    150150    INTEGER, SAVE :: ifrst = 0
    151     !$omp THREADPRIVATE(ifrst)
    152 
    153     ! =====================================================================
    154     ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
    155     ! =====================================================================
    156     ! -- Initialize mixing PDF coefficients
     151!$omp THREADPRIVATE(ifrst)
     152
     153! =====================================================================
     154! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
     155! =====================================================================
     156! -- Initialize mixing PDF coefficients
    157157    IF (ifrst == 0) THEN
    158158      ifrst = 1
     
    181181    Sigij(1:ncum, 1:nd, 1:nd) = 0.0
    182182
    183     ! =====================================================================
    184     ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
    185     ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    186     ! --- FRACTION (Sij)
    187     ! =====================================================================
     183! =====================================================================
     184! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
     185! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
     186! --- FRACTION (Sij)
     187! =====================================================================
    188188
    189189    DO i = minorig + 1, nl
     
    253253      ENDIF
    254254
    255       ! ***   if no air can entrain at level i assume that updraft detrains  ***
    256       ! ***   at that level and calculate detrained air flux and properties  ***
     255! ***   if no air can entrain at level i assume that updraft detrains  ***
     256! ***   at that level and calculate detrained air flux and properties  ***
    257257
    258258      DO il = 1, ncum
     
    282282    END DO
    283283
    284     ! =====================================================================
    285     ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
    286     ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
    287     ! =====================================================================
    288 
    289     CALL zilch(csum, nloc*nd)
     284! =====================================================================
     285! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
     286! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
     287! =====================================================================
    290288
    291289    DO il = 1, ncum
     
    293291    END DO
    294292
    295     ! ---------------------------------------------------------------
     293! ---------------------------------------------------------------
    296294    DO i = minorig + 1, nl      !Loop on origin level "i"
    297       ! ---------------------------------------------------------------
     295! ---------------------------------------------------------------
    298296
    299297      DO il = 1, ncum
     
    341339            Sbef(il) = max(0., signhpmh(il))
    342340          endif
    343         END IF
    344       END DO
    345 
    346       DO il = 1, ncum
    347         DO j = minorig, nl !Loop on destination level "j"
    348           IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    349               j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    350               lwork(il) .AND. Sij(il, i, j) > 0.0) THEN
    351             ! -----------------------------------------------
    352             IF (j > i) THEN
    353               Smid = min(Sij(il, i, j), Scrit(il))
    354               Sjmax(il) = Smid
    355               Sjmin(il) = Smid
    356               IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN
    357                 smin(il) = min(smin(il), Smid)
    358                 Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il))
    359                 Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
    360                 Sjmin(il) = min(Sjmin(il), Scrit(il))
    361                 Sbef(il) = Sij(il, i, j)
     341
     342          DO j = (icb(il) - 1), inb(il) !Loop on destination level "j"
     343            IF (lwork(il) .AND. Sij(il, i, j) > 0.0) THEN
     344              ! -----------------------------------------------
     345              IF (j > i) THEN
     346                Smid = min(Sij(il, i, j), Scrit(il))
     347                Sjmax(il) = Smid
     348                Sjmin(il) = Smid
     349                IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN
     350                  smin(il) = min(smin(il), Smid)
     351                  Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il))
     352                  Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
     353                  Sjmin(il) = min(Sjmin(il), Scrit(il))
     354                  Sbef(il) = Sij(il, i, j)
     355                END IF
     356              ELSE IF (j == i) THEN
     357                Smid = 1.
     358                Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + &
     359                            min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il))
     360                Sjmin(il) = max(Sjmin(il), sup(il))
     361                Sjmax(il) = 1.
     362                ! preparation des variables Scrit, Smin et Sbef pour la partie j>i
     363                Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
     364
     365                smin(il) = 1.
     366                Sbef(il) = max(0., signhpmh(il))
     367                supmax(il, i) = sign(Scrit(il), -signhpmh(il))
     368              ELSE IF (j < i) THEN
     369                Smid = max(Sij(il, i, j), Scrit(il))
     370                Sjmax(il) = Smid
     371                Sjmin(il) = Smid
     372                IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN
     373                  smax(il) = Smid
     374                  Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
     375                  Sjmax(il) = max(Sjmax(il), Scrit(il))
     376                  Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
     377                  Sjmin(il) = max(Sjmin(il), Scrit(il))
     378                  Sbef(il) = Sij(il, i, j)
     379                END IF
     380                IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
     381                  sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
    362382              END IF
    363             ELSE IF (j == i) THEN
    364               Smid = 1.
    365               Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + &
    366                           min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il))
    367               Sjmin(il) = max(Sjmin(il), sup(il))
    368               Sjmax(il) = 1.
    369               ! preparation des variables Scrit, Smin et Sbef pour la partie j>i
    370               Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
    371 
    372               smin(il) = 1.
    373               Sbef(il) = max(0., signhpmh(il))
    374               supmax(il, i) = sign(Scrit(il), -signhpmh(il))
    375             ELSE IF (j < i) THEN
    376               Smid = max(Sij(il, i, j), Scrit(il))
    377               Sjmax(il) = Smid
    378               Sjmin(il) = Smid
    379               IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN
    380                 smax(il) = Smid
    381                 Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
    382                 Sjmax(il) = max(Sjmax(il), Scrit(il))
    383                 Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
    384                 Sjmin(il) = max(Sjmin(il), Scrit(il))
    385                 Sbef(il) = Sij(il, i, j)
     383              ! -----------------------------------------------
     384
     385              rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
     386              Qmixmax(il) = Qmix(Sjmax(il))
     387              Qmixmin(il) = Qmix(Sjmin(il))
     388              Rmixmax(il) = Rmix(Sjmax(il))
     389              Rmixmin(il) = Rmix(Sjmin(il))
     390              sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
     391              sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
     392
     393              Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j)
     394              IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN
     395                Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il))
     396              ELSE
     397                Sigij(il, i, j) = 0.
    386398              END IF
    387               IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
    388                 sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
     399
     400              ! Compute Qent, uent, vent according to the true mixing fraction
     401              Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i)
     402              uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i)
     403              vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i)
     404
     405              ! Compute liquid water static energy of mixed draughts
     406              hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
     407              !  Heat capacity of mixed draught
     408              cpm = cpd + Qent(il, i, j)*(cpv - cpd)
     409              IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
     410                elij(il, i, j) = Qent(il, i, j) - rs(il, j)
     411                elij(il, i, j) = elij(il, i, j) + &
     412                                 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
     413                                 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
     414                elij(il, i, j) = elij(il, i, j)/ &
     415                                 (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ &
     416                                  (cpm*rrv*t(il, j)*t(il, j)))
     417              ELSE
     418                elij(il, i, j) = Qent(il, i, j) - rs(il, j)
     419                elij(il, i, j) = elij(il, i, j) + &
     420                                 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
     421                                 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
     422                elij(il, i, j) = elij(il, i, j)/ &
     423                                 (1.+lv(il, j)*lv(il, j)*rs(il, j)/ &
     424                                  (cpm*rrv*t(il, j)*t(il, j)))
     425              ENDIF
     426              elij(il, i, j) = max(elij(il, i, j), 0.)
     427
     428              elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
     429
     430              IF (j > i) THEN
     431                awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
     432                awat = amax1(awat, 0.0)
     433              ELSE
     434                awat = 0.
     435              END IF
     436
     437              ! Mixed draught temperature at level j
     438              Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j))
     439              IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
     440                hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat
     441              ELSE
     442                hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat
     443              ENDIF
     444
     445              ! ASij is the integral of P(F) over the relevant F interval
     446              ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
     447                                        Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
     448
     449              ! If I=J (detrainement and entrainement at the same level), then only the
     450              ! adiabatic ascent part of the mixture is considered
     451              IF (i == j) THEN
     452                rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
     453                Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
     454                                     Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
     455                Qent(il, i, i) = rti
     456                uent(il, i, i) = unk(il)
     457                vent(il, i, i) = vnk(il)
     458                hent(il, i, i) = hp(il, i)
     459                elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
     460                Sigij(il, i, i) = 0.
     461              END IF
    389462            END IF
    390             ! -----------------------------------------------
    391 
    392             rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
    393             Qmixmax(il) = Qmix(Sjmax(il))
    394             Qmixmin(il) = Qmix(Sjmin(il))
    395             Rmixmax(il) = Rmix(Sjmax(il))
    396             Rmixmin(il) = Rmix(Sjmin(il))
    397             sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
    398             sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
    399 
    400             Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j)
    401             IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN
    402               Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il))
    403             ELSE
    404               Sigij(il, i, j) = 0.
    405             END IF
    406 
    407             ! Compute Qent, uent, vent according to the true mixing fraction
    408             Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i)
    409             uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i)
    410             vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i)
    411 
    412             ! Compute liquid water static energy of mixed draughts
    413             hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
    414             !  Heat capacity of mixed draught
    415             cpm = cpd + Qent(il, i, j)*(cpv - cpd)
    416             IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
    417               elij(il, i, j) = Qent(il, i, j) - rs(il, j)
    418               elij(il, i, j) = elij(il, i, j) + &
    419                                (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
    420                                rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
    421               elij(il, i, j) = elij(il, i, j)/ &
    422                                (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ &
    423                                 (cpm*rrv*t(il, j)*t(il, j)))
    424             ELSE
    425               elij(il, i, j) = Qent(il, i, j) - rs(il, j)
    426               elij(il, i, j) = elij(il, i, j) + &
    427                                (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
    428                                rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
    429               elij(il, i, j) = elij(il, i, j)/ &
    430                                (1.+lv(il, j)*lv(il, j)*rs(il, j)/ &
    431                                 (cpm*rrv*t(il, j)*t(il, j)))
    432             ENDIF
    433             elij(il, i, j) = max(elij(il, i, j), 0.)
    434 
    435             elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
    436 
    437             IF (j > i) THEN
    438               awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
    439               awat = amax1(awat, 0.0)
    440             ELSE
    441               awat = 0.
    442             END IF
    443 
    444             ! Mixed draught temperature at level j
    445             IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
    446               Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j))
    447               hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat
    448             ELSE
    449               Tm = t(il, j) + (Qent(il, i, j) - elij(il, i, j) - rs(il, j))*rrv*t(il, j)*t(il, j)/(lv(il, j)*rs(il, j))
    450               hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat
    451             ENDIF
    452 
    453             ! ASij is the integral of P(F) over the relevant F interval
    454             ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
    455                                       Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
    456 
    457             ! If I=J (detrainement and entrainement at the same level), then only the
    458             ! adiabatic ascent part of the mixture is considered
    459             IF (i == j) THEN
    460               rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
    461               Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
    462                                    Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
    463               Qent(il, i, i) = rti
     463          END DO ! Loop j = (icb(il) - 1), inb(il)
     464
     465          IF (lwork(il)) THEN
     466            csum(il, i) = 0
     467            DO j = (icb(il) - 1), inb(il)
     468              ASij(il) = amax1(1.0E-16, ASij(il))
     469              ASij_inv(il) = 1.0/ASij(il)
     470              ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
     471              IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
     472              Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
     473              csum(il, i) = csum(il, i) + Ment(il, i, j)
     474            END DO
     475            IF (csum(il, i) < 1.) THEN
     476              nent(il, i) = 0
     477              Ment(il, i, i) = 1.
     478              Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
    464479              uent(il, i, i) = unk(il)
    465480              vent(il, i, i) = vnk(il)
    466               hent(il, i, i) = hp(il, i)
    467481              elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
    468               Sigij(il, i, i) = 0.
    469             END IF
    470           END IF
    471         END DO ! End loop on destination level "j"
    472         ! ---------------------------------------------------------------
    473       END DO
    474 
    475       DO il = 1, ncum
    476         IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN
    477           ASij(il) = amax1(1.0E-16, ASij(il))
    478           ASij_inv(il) = 1.0/ASij(il)
    479           !   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
    480           IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
    481           csum(il, i) = 0.0
    482         END IF
    483       END DO
    484 
    485       DO j = minorig, nl
    486         DO il = 1, ncum
    487           IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
    488               j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
    489             Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
    490           END IF
    491         END DO
    492       END DO
    493 
    494       DO j = minorig, nl
    495         DO il = 1, ncum
    496           IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
    497               j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
    498             csum(il, i) = csum(il, i) + Ment(il, i, j)
    499           END IF
    500         END DO
    501       END DO
    502 
    503       DO il = 1, ncum
    504         IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN
    505           nent(il, i) = 0
    506           Ment(il, i, i) = 1.
    507           Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
    508           uent(il, i, i) = unk(il)
    509           vent(il, i, i) = vnk(il)
    510           elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
    511           IF (fl_cor_ebil .GE. 2) THEN
    512             hent(il, i, i) = hp(il, i)
    513             Sigij(il, i, i) = 0.0
    514           ELSE
    515             Sij(il, i, i) = 0.0
    516           ENDIF
    517         END IF
    518       END DO ! il
    519 
    520       ! ---------------------------------------------------------------
    521     END DO              ! End loop on origin level "i"
    522     ! ---------------------------------------------------------------
    523     RETURN
     482              IF (fl_cor_ebil .GE. 2) THEN
     483                hent(il, i, i) = hp(il, i)
     484                Sigij(il, i, i) = 0.0
     485              ELSE
     486                Sij(il, i, i) = 0.0
     487              ENDIF
     488            ENDIF
     489          END IF
     490        END IF ! i >= icb(il) .AND. i <= inb(il)
     491      END DO ! Loop il = 1, ncum
     492    END DO ! End loop on origin level "i"
    524493  END SUBROUTINE
    525494
Note: See TracChangeset for help on using the changeset viewer.