Changeset 3710 for LMDZ6/branches


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

cv3p_mixing_new : swapped loops ij and j

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90

    r3709 r3710  
    106106    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: hent
    107107    INTEGER, DIMENSION(nloc, nd), INTENT(OUT)        :: nent
    108  
    109   !local variables:
     108
     109!local variables:
    110110    INTEGER i, j, k, il, im, jm
    111111    REAL Smid
    … …  
    126126    REAL                               :: Tm         !Mixed draught temperature
    127127    LOGICAL, DIMENSION(nloc)          :: lwork
    128  
     128
    129129    REAL amxupcrit, df, ff
    130130    INTEGER nstep
    131  
     131
    132132    INTEGER, SAVE                                       :: igout = 1
    133   !$omp THREADPRIVATE(igout)
    134  
    135   ! --   Mixing probability distribution functions
    136  
     133    !$omp THREADPRIVATE(igout)
     134
     135    ! --   Mixing probability distribution functions
     136
    137137    REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
    138  
     138
    139139    Qcoef1(F) = tanh(F/gammas)
    140140    Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
    … …  
    147147    Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
    148148    Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
    149  
     149
    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
    … …  
    162162        fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
    163163    END IF
    164  
     164
    165165    nent(:ncum, :nl) = 0
    166166    elij(:ncum, :nl, :nl) = 0
    167167    hent(:ncum, :nl, :nl) = 0
    168  
     168
    169169    DO j = 1, nl
    170170      DO k = 1, nl
    … …  
    176176      END DO
    177177    END DO
    178  
     178
    179179    Ment(1:ncum, 1:nd, 1:nd) = 0.0
    180180    Sij(1:ncum, 1:nd, 1:nd) = 0.0
    181181    Sigij(1:ncum, 1:nd, 1:nd) = 0.0
    182  
    183   ! =====================================================================
    184   ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
    185   ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    186   ! --- FRACTION (Sij)
    187   ! =====================================================================
    188  
     182
     183    ! =====================================================================
     184    ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
     185    ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
     186    ! --- FRACTION (Sij)
     187    ! =====================================================================
     188
    189189    DO i = minorig + 1, nl
    190190      IF (ok_entrain) THEN
    … …  
    230230                nent(il, i) = nent(il, i) + 1
    231231              END IF
    232  
     232
    233233              Sij(il, i, j) = amax1(0.0, Sij(il, i, j))
    234234              Sij(il, i, j) = amin1(1.0, Sij(il, i, j))
    … …  
    245245        ENDDO
    246246      ENDIF ! (ok_entrain)
    247  
     247
    248248      IF (prt_level >= 10) THEN
    249249        print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)
    … …  
    252252        ENDIF
    253253      ENDIF
    254  
    255   ! ***   if no air can entrain at level i assume that updraft detrains  ***
    256   ! ***   at that level and calculate detrained air flux and properties  ***
    257  
     254
     255      ! ***   if no air can entrain at level i assume that updraft detrains  ***
     256      ! ***   at that level and calculate detrained air flux and properties  ***
     257
    258258      DO il = 1, ncum
    259259        IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN
    … …  
    270270      END DO
    271271    END DO ! i = minorig + 1, nl
    272  
     272
    273273    DO j = minorig, nl
    274274      DO i = minorig, nl
    … …  
    281281      END DO
    282282    END DO
    283  
    284   ! =====================================================================
    285   ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
    286   ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
    287   ! =====================================================================
    288  
     283
     284    ! =====================================================================
     285    ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
     286    ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
     287    ! =====================================================================
     288
    289289    CALL zilch(csum, nloc*nd)
    290  
     290
    291291    DO il = 1, ncum
    292292      lwork(il) = .FALSE.
    293293    END DO
    294  
    295   ! ---------------------------------------------------------------
     294
     295    ! ---------------------------------------------------------------
    296296    DO i = minorig + 1, nl      !Loop on origin level "i"
    297   ! ---------------------------------------------------------------
    298      
     297      ! ---------------------------------------------------------------
     298
    299299      DO il = 1, ncum
    300300        IF (i >= icb(il) .AND. i <= inb(il)) THEN
    301301          signhpmh(il) = sign(1., hp(il, i) - h(il, i))
    302          
    303           Sx(il) = max( maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il) )
    304        
     302
     303          Sx(il) = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh(il))
     304
    305305          lwork(il) = (nent(il, i) /= 0)
    306306          rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
    … …  
    319319          Scrit(il) = min(anum/denom, 1.)
    320320          alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti)
    321  
     321
    322322          ! Find new critical value Scrit2
    323323          ! such that : Sij > Scrit2  => mixed draught will detrain at J<I
    … …  
    326326                   Scrit(il)*max(0., signhpmh(il))
    327327          Scrit(il) = Scrit2
    328  
     328
    329329          ! Correction pour la nouvelle logique; la correction pour ALT
    330330          ! est un peu au hazard
    331331          IF (Scrit(il) <= 0.0) Scrit(il) = 0.0
    332332          IF (alt <= 0.0) Scrit(il) = 1.0
    333  
     333
    334334          smax(il) = 0.0
    335335          ASij(il) = 0.0
    336336          sup(il) = 0.      ! upper S-value reached by descending draughts
    337          
    338           if( i < inb(il)) then
     337
     338          if (i < inb(il)) then
    339339            Sbef(il) = Sij(il, i, inb(il))
    340340          else
    … …  
    343343        END IF
    344344      END DO
    345  
    346   ! ---------------------------------------------------------------
    347       DO j = minorig, nl         !Loop on destination level "j"
    348   ! ---------------------------------------------------------------       
    349   ! -----------------------------------------------
    350         IF (j > i) THEN
    351   ! -----------------------------------------------
    352           DO il = 1, ncum
    353             IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    354                 j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    355                 lwork(il)) THEN
    356               IF (Sij(il, i, j) > 0.0) THEN
    357                 Smid = min(Sij(il, i, j), Scrit(il))
    358                 Sjmax(il) = Smid
    359                 Sjmin(il) = Smid
    360                 IF (Smid < smin(il) .AND. Sij(il, i, j + 1) < Smid) THEN
    361                   smin(il) = min( smin(il), Smid )
    362                   Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il))
    363                   Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
    364                   Sjmin(il) = min(Sjmin(il), Scrit(il))
    365                   Sbef(il) = Sij(il, i, j)
    366                 END IF
    367               END IF
    368             END IF
    369           END DO
    370   ! -----------------------------------------------
    371         ELSE IF (j == i) THEN
    372   ! -----------------------------------------------
    373           DO il = 1, ncum
    374             IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    375                 j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    376                 lwork(il)) THEN
    377               IF (Sij(il, i, j) > 0.0) THEN
    378                 Smid = 1.
    379                 Sjmin(il) = max((Sij(il, i, j - 1) + Smid)/2., Scrit(il))*max(0., -signhpmh(il)) + &
    380                             min((Sij(il, i, j + 1) + Smid)/2., Scrit(il))*max(0., signhpmh(il))
    381                 Sjmin(il) = max(Sjmin(il), sup(il))
    382                 Sjmax(il) = 1.
    383   ! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
    384                 Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
    385  
    386                 smin(il) = 1.
    387                 Sbef(il) = max(0., signhpmh(il))
    388                 supmax(il, i) = sign(Scrit(il), -signhpmh(il))
    389               END IF
    390             END IF
    391           END DO
    392   ! -----------------------------------------------
    393         ELSE IF (j < i) THEN
    394   ! -----------------------------------------------
    395           DO il = 1, ncum
    396             IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    397                 j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    398                 lwork(il)) THEN
    399               IF (Sij(il, i, j) > 0.0) THEN
    400                 Smid = max(Sij(il, i, j), Scrit(il))
    401                 Sjmax(il) = Smid
    402                 Sjmin(il) = Smid
    403                 IF (Smid > smax(il) .AND. Sij(il, i, j + 1) > Smid) THEN
    404                   smax(il) = Smid
    405                   Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
    406                   Sjmax(il) = max(Sjmax(il), Scrit(il))
    407                   Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
    408                   Sjmin(il) = max(Sjmin(il), Scrit(il))
    409                   Sbef(il) = Sij(il, i, j)
    410                 END IF
    411                 IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
    412                   sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
    413               END IF
    414             END IF
    415           END DO
    416   ! -----------------------------------------------
    417         END IF
    418   ! -----------------------------------------------
    419  
    420         DO il = 1, ncum
     345
     346      DO il = 1, ncum
     347        DO j = minorig, nl !Loop on destination level "j"
    421348          IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    422349              j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    423               lwork(il)) THEN
    424             IF (Sij(il, i, j) > 0.0) THEN
     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)
     362              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)
     386              END IF
     387              IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
     388                sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
     389            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
    425460              rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
    426               Qmixmax(il) = Qmix(Sjmax(il))
    427               Qmixmin(il) = Qmix(Sjmin(il))
    428               Rmixmax(il) = Rmix(Sjmax(il))
    429               Rmixmin(il) = Rmix(Sjmin(il))
    430               sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
    431               sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
    432  
    433               Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j)
    434               IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN
    435                 Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il))
    436               ELSE
    437                 Sigij(il, i, j) = 0.
    438               END IF
    439  
    440   ! --    Compute Qent, uent, vent according to the true mixing fraction
    441               Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i)
    442               uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i)
    443               vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i)
    444  
    445   ! --     Compute liquid water static energy of mixed draughts
    446               hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
    447   !  Heat capacity of mixed draught
    448               cpm = cpd + Qent(il, i, j)*(cpv - cpd)
    449               IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
    450                 elij(il, i, j) = Qent(il, i, j) - rs(il, j)
    451                 elij(il, i, j) = elij(il, i, j) + &
    452                                  (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
    453                                  rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
    454                 elij(il, i, j) = elij(il, i, j)/ &
    455                                  (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ &
    456                                   (cpm*rrv*t(il, j)*t(il, j)))
    457               ELSE
    458                 elij(il, i, j) = Qent(il, i, j) - rs(il, j)
    459                 elij(il, i, j) = elij(il, i, j) + &
    460                                  (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
    461                                  rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
    462                 elij(il, i, j) = elij(il, i, j)/ &
    463                                  (1.+lv(il, j)*lv(il, j)*rs(il, j)/ &
    464                                   (cpm*rrv*t(il, j)*t(il, j)))
    465               ENDIF
    466               elij(il, i, j) = max(elij(il, i, j), 0.)
    467  
    468               elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
    469  
    470               IF (j > i) THEN
    471                 awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
    472                 awat = amax1(awat, 0.0)
    473               ELSE
    474                 awat = 0.
    475               END IF
    476  
    477   ! Mixed draught temperature at level j
    478               IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
    479                 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))
    480                 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat
    481               ELSE
    482                 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))
    483                 hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat
    484               ENDIF
    485  
    486   ! --      ASij is the integral of P(F) over the relevant F interval
    487               ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
    488                                         Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
    489  
     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
     464              uent(il, i, i) = unk(il)
     465              vent(il, i, i) = vnk(il)
     466              hent(il, i, i) = hp(il, i)
     467              elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
     468              Sigij(il, i, i) = 0.
    490469            END IF
    491470          END IF
    492         END DO
    493  
    494   ! --    If I=J (detrainement and entrainement at the same level), then only the
    495   ! --    adiabatic ascent part of the mixture is considered
    496         IF (i == j) THEN
    497           DO il = 1, ncum
    498             IF (i >= icb(il) .AND. i <= inb(il) .AND. &
    499                 j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
    500                 lwork(il)) THEN
    501               IF (Sij(il, i, j) > 0.0) THEN
    502                 rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
    503                 Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
    504                                      Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
    505                 Qent(il, i, i) = rti
    506                 uent(il, i, i) = unk(il)
    507                 vent(il, i, i) = vnk(il)
    508                 hent(il, i, i) = hp(il, i)
    509                 elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
    510                 Sigij(il, i, i) = 0.
    511               END IF
    512             END IF
    513           END DO
    514         END IF
    515  
    516   ! ---------------------------------------------------------------
    517   175 END DO        ! End loop on destination level "j"
    518   ! ---------------------------------------------------------------
    519  
     471        END DO ! End loop on destination level "j"
     472        ! ---------------------------------------------------------------
     473      END DO
     474
    520475      DO il = 1, ncum
    521476        IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN
    522477          ASij(il) = amax1(1.0E-16, ASij(il))
    523478          ASij_inv(il) = 1.0/ASij(il)
    524   !   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
     479          !   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
    525480          IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
    526481          csum(il, i) = 0.0
    527482        END IF
    528483      END DO
    529  
     484
    530485      DO j = minorig, nl
    531486        DO il = 1, ncum
    … …  
    536491        END DO
    537492      END DO
    538  
     493
    539494      DO j = minorig, nl
    540495        DO il = 1, ncum
    … …  
    545500        END DO
    546501      END DO
    547  
     502
    548503      DO il = 1, ncum
    549504        IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN
    … …  
    562517        END IF
    563518      END DO ! il
    564  
    565   ! ---------------------------------------------------------------
    566   789 END DO              ! End loop on origin level "i"
    567   ! ---------------------------------------------------------------
     519
     520      ! ---------------------------------------------------------------
     521    END DO              ! End loop on origin level "i"
     522    ! ---------------------------------------------------------------
    568523    RETURN
    569524  END SUBROUTINE
Note: See TracChangeset for help on using the changeset viewer.