source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/cv3p_mixing.f90 @ 4703

Last change on this file since 4703 was 3717, checked in by adurocher, 5 years ago

cv3p_mixing_new : VLAs removed from block construct (intel 19 bug)

Intel compiler 19 crashes (internal compiler error) when VLA (variable length arrays) are used within block constructs
File size: 47.4 KB
Line 
1module cv3p_mixing_mod
2  implicit none
3
4contains
5
6  ! **************************************************************
7  ! *
8  ! CV3P_MIXING : compute mixed draught properties and,         *
9  ! within a scaling factor, mixed draught        *
10  ! mass fluxes.                                  *
11  ! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
12  ! modified by :                                               *
13  ! **************************************************************
14  SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
15                         ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
16                         unk, vnk, hp, tv, tvp, ep, clw, sig, &
17                         Ment, Qent, hent, uent, vent, nent, &
18                         Sigij, elij, supmax, Ments, Qents, traent)
19
20    !inputs:
21    INTEGER, INTENT(IN)                               :: ncum, nd, na
22    INTEGER, INTENT(IN)                               :: ntra, nloc
23    INTEGER, DIMENSION(nloc), INTENT(IN)             :: icb, inb, nk
24    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: sig
25    REAL, DIMENSION(nloc), INTENT(IN)                :: unk, vnk
26    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: qta
27    REAL, DIMENSION(nloc, nd + 1), INTENT(IN)          :: ph
28    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: t, rr, rs
29    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: u, v
30    REAL, DIMENSION(nloc, nd, ntra), INTENT(IN)      :: tra  ! input of convect3
31    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lv
32    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lf
33    REAL, DIMENSION(nloc, na), INTENT(IN)            :: frac  !ice fraction in condensate
34    REAL, DIMENSION(nloc, na), INTENT(IN)            :: h   !liquid water static energy of environment
35    REAL, DIMENSION(nloc, na), INTENT(IN)            :: hp  !liquid water static energy of air shed from adiab. asc.
36    REAL, DIMENSION(nloc, na), INTENT(IN)            :: tv, tvp
37    REAL, DIMENSION(nloc, na), INTENT(IN)            :: ep, clw
38
39    !outputs:
40    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Ment, Qent
41    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: uent, vent
42    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Sigij, elij
43    REAL, DIMENSION(nloc, na), INTENT(OUT)           :: supmax            ! Highest mixing fraction of mixed
44    ! updraughts with the sign of (h-hp)
45    REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent
46    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: Ments, Qents
47    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: hent
48    INTEGER, DIMENSION(nloc, nd), INTENT(OUT)        :: nent
49
50    logical, parameter :: use_old = .false.
51
52    if (use_old) then
53      call cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, &
54                           ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
55                           unk, vnk, hp, tv, tvp, ep, clw, sig, &
56                           Ment, Qent, hent, uent, vent, nent, &
57                           Sigij, elij, supmax, Ments, Qents, traent)
58    else
59      call cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, &
60                           t, rr, rs, u, v, h, lv, lf, frac, qta, &
61                           unk, vnk, hp, ep, clw, &
62                           Ment, Qent, hent, uent, vent, nent, &
63                           Sigij, elij, supmax)
64    endif
65
66  END SUBROUTINE
67
68  ! modified by : Arnaud Durocher (04/2020) : changed loops order
69  SUBROUTINE cv3p_mixing_new(nloc, ncum, nd, na, icb, inb, &
70                             t, rr, rs, u, v, h, lv, lf, frac, qta, &
71                             unk, vnk, hp, ep, clw, &
72                             Ment, Qent, hent, uent, vent, nent, &
73                             Sigij, elij, supmax)
74
75    USE print_control_mod, ONLY: prt_level
76    USE add_phys_tend_mod, ONLY: fl_cor_ebil
77
78    IMPLICIT NONE
79
80    include "cvthermo.h"
81    include "cv3param.h"
82    include "YOMCST2.h"
83    include "cvflag.h"
84
85    !inputs:
86    INTEGER, INTENT(IN)                              :: ncum, nd, na
87    INTEGER, INTENT(IN)                              :: nloc
88    INTEGER, DIMENSION(nloc), INTENT(IN)             :: icb, inb
89    REAL, DIMENSION(nloc), INTENT(IN)                :: unk, vnk
90    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: qta
91    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: t, rr, rs
92    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: u, v
93    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lv
94    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lf
95    REAL, DIMENSION(nloc, na), INTENT(IN)            :: frac  !ice fraction in condensate
96    REAL, DIMENSION(nloc, na), INTENT(IN)            :: h   !liquid water static energy of environment
97    REAL, DIMENSION(nloc, na), INTENT(IN)            :: hp  !liquid water static energy of air shed from adiab. asc.
98    REAL, DIMENSION(nloc, na), INTENT(IN)            :: ep, clw
99
100    !outputs:
101    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Ment, Qent
102    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: uent, vent
103    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Sigij, elij
104    REAL, DIMENSION(nloc, na), INTENT(OUT)           :: supmax            ! Highest mixing fraction of mixed
105    ! updraughts with the sign of (h-hp)
106    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: hent
107    INTEGER, DIMENSION(nloc, nd), INTENT(OUT)        :: nent
108
109!local variables:
110    INTEGER i, j, k, il
111    REAL, DIMENSION(nloc, nd)          :: ASij
112    REAL, DIMENSION(nloc, nd, nd)      :: Sij
113    real :: Sjmin(nd), Sjmax(nd), Qmixmax(nd), Qmixmin(nd), Rmixmax(nd), Rmixmin(nd)
114
115    INTEGER, SAVE                      :: igout = 1
116    !$omp THREADPRIVATE(igout)
117
118    ! --   Mixing probability distribution functions
119
120    REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
121
122    Qcoef1(F) = tanh(F/gammas)
123    Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
124    QFF(F) = max(min(F, 1.), 0.)
125    QFFf(F) = min(QFF(F), scut)
126    Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max
127    Rmix1(F) = (gammas*log(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)*Qcoef1max)/Qcoef2max
128    Qmix2(F) = -log(1.-QFFf(F))/scut
129    Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut
130    Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
131    Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
132
133    INTEGER, SAVE :: ifrst = 0
134    !$omp THREADPRIVATE(ifrst)
135
136    ! =====================================================================
137    ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
138    ! =====================================================================
139    ! -- Initialize mixing PDF coefficients
140    IF (ifrst == 0) THEN
141      ifrst = 1
142      Qcoef1max = Qcoef1(Fmax)
143      Qcoef2max = Qcoef2(Fmax)
144      print *, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
145        fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
146    END IF
147
148    nent(:ncum, :nl) = 0
149    elij(:ncum, :nl, :nl) = 0
150    hent(:ncum, :nl, :nl) = 0
151
152    DO j = 1, nl
153      DO k = 1, nl
154        DO i = 1, ncum
155          Qent(i, k, j) = rr(i, j)
156          uent(i, k, j) = u(i, j)
157          vent(i, k, j) = v(i, j)
158        END DO
159      END DO
160    END DO
161
162    Ment(1:ncum, 1:nd, 1:nd) = 0.0
163    Sij(1:ncum, 1:nd, 1:nd) = 0.0
164
165    ! =====================================================================
166    ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
167    ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
168    ! --- FRACTION (Sij)
169    ! =====================================================================
170
171    IF (ok_entrain) THEN
172      ! NOTE : (04/2020 AD) this loop order gave good results in cv3a_tracer
173      !        if( icb(il) <= i <= inb(il) ) is probably more expensive than reading non-contiguously
174      DO il = 1, ncum
175        DO i = icb(il), inb(il)
176          DO j = (icb(il) - 1), inb(il)
177            block
178              real :: rti, bf2, anum, denom, dei, altem, cwat, stemp
179
180              rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
181              bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd)
182              IF (cvflag_ice .AND. t(il, j) <= 263.15) THEN
183                bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* &
184                                                  lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd)
185              END IF
186              anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j))
187            denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j)
188            dei = denom
189            IF (abs(dei) < 0.01) dei = 0.01
190            Sij(il, i, j) = anum/dei
191            Sij(il, i, i) = 1.0
192            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j)
193            altem = altem/bf2
194            cwat = clw(il, j)*(1.-ep(il, j))
195              stemp = Sij(il, i, j)
196              IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN
197                IF (cvflag_ice) THEN
198                  anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2)
199                  denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti)
200                ELSE
201                  anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2)
202                  denom = denom + lv(il, j)*(rr(il, i) - rti)
203                END IF
204                IF (abs(denom) < 0.01) denom = 0.01
205                Sij(il, i, j) = anum/denom
206                altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j)
207                altem = altem - (bf2 - 1.)*cwat
208              END IF
209              IF (Sij(il, i, j) > 0.0) THEN
210                Ment(il, i, j) = 1.
211                elij(il, i, j) = altem
212                elij(il, i, j) = amax1(0.0, elij(il, i, j))
213                nent(il, i) = nent(il, i) + 1
214              END IF
215
216              Sij(il, i, j) = amax1(0.0, Sij(il, i, j))
217              Sij(il, i, j) = amin1(1.0, Sij(il, i, j))
218            end block
219          END DO
220          IF (nent(il, i) == 0) THEN
221            Ment(il, i, i) = 1.
222            Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
223            uent(il, i, i) = unk(il)
224            vent(il, i, i) = vnk(il)
225            IF (fl_cor_ebil .GE. 2) THEN
226              hent(il, i, i) = hp(il, i)
227            ENDIF
228            elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
229            Sij(il, i, i) = 0.0
230          END IF
231        END DO
232      END DO
233    ELSE
234      DO i = minorig + 1, nl
235        DO il = 1, ncum
236          IF ((i >= icb(il)) .AND. (i <= inb(il))) THEN
237            Ment(il, i, i) = 1.
238            Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
239            uent(il, i, i) = unk(il)
240            vent(il, i, i) = vnk(il)
241            IF (fl_cor_ebil .GE. 2) THEN
242              hent(il, i, i) = hp(il, i)
243            ENDIF
244            elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
245            Sij(il, i, i) = 0.0
246          END IF
247        END DO
248      END DO
249    ENDIF
250
251    DO i = minorig + 1, nl
252      DO j = minorig, nl
253        IF (prt_level >= 10) THEN
254          print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)
255          !IF (nent(igout, i) .gt. 0) THEN
256          !  print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout))
257          !ENDIF
258        ENDIF
259      END DO
260    END DO
261
262    Sigij(1:ncum, 1:nd, 1:nd) = 0.0
263    DO il = 1, ncum
264      DO i = icb(il), inb(il)
265        DO j = (icb(il) - 1), inb(il)
266          Sigij(il, i, j) = Sij(il, i, j)
267        END DO
268      END DO
269    END DO
270
271    ! =====================================================================
272    ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
273    ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
274    ! =====================================================================
275    DO il = 1, ncum
276      DO i = icb(il), inb(il)      !Loop on origin level "i"
277        block
278          real :: signhpmh, Scrit
279
280          block
281            real :: rti, anum, denom, alt, Sx
282
283            signhpmh = sign(1., hp(il, i) - h(il, i))
284
285            rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
286            IF (cvflag_ice) THEN
287              anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* &
288                     (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i))
289              denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* &
290                      (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti)
291            ELSE
292              anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + &
293                     (cpv - cpd)*t(il, i)*(rti - rr(il, i))
294              denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + &
295                      (cpd - cpv)*t(il, i)*(rr(il, i) - rti)
296            END IF
297            IF (abs(denom) < 0.01) denom = 0.01
298            Scrit = min(anum/denom, 1.)
299            alt = rti - rs(il, i) + Scrit*(rr(il, i) - rti)
300
301            if (alt <= 0) then
302              Scrit = 1.0
303            else
304              ! Get max of Sij
305              Sx = max(maxval(Sij(il, i, i + 1:inb(il))), 0., signhpmh)
306              ! Find new critical value Scrit
307              ! such that : Sij > Scrit  => mixed draught will detrain at J<I
308              !             Sij < Scrit  => mixed draught will detrain at J>I
309              Scrit = max(0., min(Scrit, Sx*max(0., -signhpmh)) + Scrit*max(0., signhpmh))
310            endif
311
312          end block
313
314          ASij(il, i) = 0.0
315
316          ! Compute all Sjmax(j) and Sjmin(j)
317          block
318            real :: smin, smax, sup, Sbef, Smid
319            smin = 1
320            smax = 0.0
321            sup = 0. ! upper S-value reached by descending draughts
322            ! Glitchy : why so complicated?
323            if (i < inb(il)) then
324              Sbef = Sij(il, i, inb(il))
325            else
326              Sbef = max(0., signhpmh)
327            endif
328
329            DO j = (icb(il) - 1), i - 1
330              IF (Sij(il, i, j) > 0.0) THEN
331                Smid = max(Sij(il, i, j), Scrit)
332                Sjmax(j) = Smid
333                Sjmin(j) = Smid
334                IF (Smid > smax .AND. Sij(il, i, j + 1) > Smid) THEN
335                  smax = Smid
336                  Sjmax(j) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
337                  Sjmax(j) = max(Sjmax(j), Scrit)
338                  Sjmin(j) = min((Sbef + Sij(il, i, j))/2., Sij(il, i, j))
339                  Sjmin(j) = max(Sjmin(j), Scrit)
340                  Sbef = Sij(il, i, j)
341                END IF
342                IF (abs(Sjmin(j) - Sjmax(j)) > 1.E-10) &
343                  sup = max(Sjmin(j), Sjmax(j), sup)
344              ENDIF
345            END DO
346            j = i
347            IF (Sij(il, i, j) > 0.0) THEN
348              Smid = 1.
349              Sjmin(j) = max((Sij(il, i, j - 1) + Smid)/2., Scrit)*max(0., -signhpmh) + &
350                         min((Sij(il, i, j + 1) + Smid)/2., Scrit)*max(0., signhpmh)
351              Sjmin(j) = max(Sjmin(j), sup)
352              Sjmax(j) = 1.
353              ! preparation des variables Scrit, Smin et Sbef pour la partie j>i
354              Scrit = min(Sjmin(j), Sjmax(j), Scrit)
355              Sbef = max(0., signhpmh)
356              supmax(il, i) = sign(Scrit, -signhpmh)
357            ENDIF
358            DO j = i + 1, inb(il)
359              IF (Sij(il, i, j) > 0.0) THEN
360                Smid = min(Sij(il, i, j), Scrit)
361                Sjmax(j) = Smid
362                Sjmin(j) = Smid
363                IF (Smid < smin .AND. Sij(il, i, j + 1) < Smid) THEN
364                  smin = min(smin, Smid)
365                  Sjmax(j) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit)
366                  Sjmin(j) = max((Sbef + Sij(il, i, j))/2., Sij(il, i, j))
367                  Sjmin(j) = min(Sjmin(j), Scrit)
368                  Sbef = Sij(il, i, j)
369                END IF
370              ENDIF
371            END DO
372          end block
373
374          DO j = (icb(il) - 1), inb(il)
375            Qmixmax(j) = Qmix(Sjmax(j))
376            Qmixmin(j) = Qmix(Sjmin(j))
377            Rmixmax(j) = Rmix(Sjmax(j))
378            Rmixmin(j) = Rmix(Sjmin(j))
379          END DO
380
381          DO j = (icb(il) - 1), inb(il) !Loop on destination level "j"
382            IF (Sij(il, i, j) > 0.0) THEN
383              block
384                real :: sqmrmax, sqmrmin, cpm, awat, Tm, rti
385
386                rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
387
388                sqmrmax = Sjmax(j)*Qmixmax(j) - Rmixmax(j)
389                sqmrmin = Sjmin(j)*Qmixmin(j) - Rmixmin(j)
390
391                Ment(il, i, j) = abs(Qmixmax(j) - Qmixmin(j))*Ment(il, i, j)
392                IF (abs(Qmixmax(j) - Qmixmin(j)) > 1.E-10) THEN
393                  Sigij(il, i, j) = (sqmrmax - sqmrmin)/(Qmixmax(j) - Qmixmin(j))
394                ELSE
395                  Sigij(il, i, j) = 0.
396                END IF
397
398                ! Compute Qent, uent, vent according to the true mixing fraction
399                Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i)
400                uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i)
401                vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i)
402
403                ! Compute liquid water static energy of mixed draughts
404                hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
405                !  Heat capacity of mixed draught
406                cpm = cpd + Qent(il, i, j)*(cpv - cpd)
407
408                elij(il, i, j) = Qent(il, i, j) - rs(il, j)
409                elij(il, i, j) = elij(il, i, j) + &
410                                 (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
411                                 rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
412                IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
413                  elij(il, i, j) = elij(il, i, j)/ &
414                                   (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ &
415                                    (cpm*rrv*t(il, j)*t(il, j)))
416                ELSE
417
418                  elij(il, i, j) = elij(il, i, j)/ &
419                                   (1.+lv(il, j)*lv(il, j)*rs(il, j)/ &
420                                    (cpm*rrv*t(il, j)*t(il, j)))
421                ENDIF
422                elij(il, i, j) = max(elij(il, i, j), 0.)
423                elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
424
425                IF (j > i) THEN
426                  awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
427                  awat = amax1(awat, 0.0)
428                ELSE
429                  awat = 0.
430                END IF
431
432                ! Mixed draught temperature at level j
433                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))
434                IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
435                  hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat
436                ELSE
437                  hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat
438                ENDIF
439
440                ! ASij is the integral of P(F) over the relevant F interval
441                ASij(il, i) = ASij(il, i) + abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - &
442                                                Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j))
443
444                ! If I=J (detrainement and entrainement at the same level), then only the
445                ! adiabatic ascent part of the mixture is considered
446                IF (i == j) THEN
447                  rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
448                  Ment(il, i, i) = abs(Qmixmax(j)*(1.-Sjmax(j)) + Rmixmax(j) - &
449                                       Qmixmin(j)*(1.-Sjmin(j)) - Rmixmin(j))
450                  Qent(il, i, i) = rti
451                  uent(il, i, i) = unk(il)
452                  vent(il, i, i) = vnk(il)
453                  hent(il, i, i) = hp(il, i)
454                  elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
455                  Sigij(il, i, i) = 0.
456                END IF
457              end block
458            END IF
459          END DO ! Loop j = (icb(il) - 1), inb(il)
460        end block
461      END DO ! Loop i = icb(il), inb(il)
462    END DO ! Loop il = 1, ncum
463
464    DO il = 1, ncum
465      DO i = icb(il), inb(il) !Loop on origin level "i"
466        block
467          real :: csum, ASij_inv
468          csum = 0
469          DO j = (icb(il) - 1), inb(il)
470            ASij(il, i) = amax1(1.0E-16, ASij(il, i))
471            ASij_inv = 1.0/ASij(il, i)
472            ! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
473            IF (ASij_inv > 100.) ASij_inv = 0.
474            Ment(il, i, j) = Ment(il, i, j)*ASij_inv
475            csum = csum + Ment(il, i, j)
476          END DO
477          IF (csum < 1.) THEN
478            nent(il, i) = 0
479            Ment(il, i, i) = 1.
480            Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
481            uent(il, i, i) = unk(il)
482            vent(il, i, i) = vnk(il)
483            elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
484            IF (fl_cor_ebil .GE. 2) THEN
485              hent(il, i, i) = hp(il, i)
486              Sigij(il, i, i) = 0.0
487            ELSE
488              Sij(il, i, i) = 0.0
489            ENDIF
490          ENDIF
491        end block
492      END DO ! Loop il = 1, ncum
493    END DO ! End loop on origin level "i"
494  END SUBROUTINE
495
496  ! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15
497  SUBROUTINE cv3p_mixing_old(nloc, ncum, nd, na, ntra, icb, nk, inb, &
498                             ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
499                             unk, vnk, hp, tv, tvp, ep, clw, sig, &
500                             Ment, Qent, hent, uent, vent, nent, &
501                             Sigij, elij, supmax, Ments, Qents, traent)
502
503    USE print_control_mod, ONLY: prt_level
504    USE add_phys_tend_mod, ONLY: fl_cor_ebil
505
506    IMPLICIT NONE
507
508    include "cvthermo.h"
509    include "cv3param.h"
510    include "YOMCST2.h"
511    include "cvflag.h"
512
513    !inputs:
514    INTEGER, INTENT(IN)                               :: ncum, nd, na
515    INTEGER, INTENT(IN)                               :: ntra, nloc
516    INTEGER, DIMENSION(nloc), INTENT(IN)             :: icb, inb, nk
517    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: sig
518    REAL, DIMENSION(nloc), INTENT(IN)                :: unk, vnk
519    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: qta
520    REAL, DIMENSION(nloc, nd + 1), INTENT(IN)          :: ph
521    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: t, rr, rs
522    REAL, DIMENSION(nloc, nd), INTENT(IN)            :: u, v
523    REAL, DIMENSION(nloc, nd, ntra), INTENT(IN)      :: tra  ! input of convect3
524    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lv
525    REAL, DIMENSION(nloc, na), INTENT(IN)            :: lf
526    REAL, DIMENSION(nloc, na), INTENT(IN)            :: frac  !ice fraction in condensate
527    REAL, DIMENSION(nloc, na), INTENT(IN)            :: h   !liquid water static energy of environment
528    REAL, DIMENSION(nloc, na), INTENT(IN)            :: hp  !liquid water static energy of air shed from adiab. asc.
529    REAL, DIMENSION(nloc, na), INTENT(IN)            :: tv, tvp
530    REAL, DIMENSION(nloc, na), INTENT(IN)            :: ep, clw
531
532    !outputs:
533    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Ment, Qent
534    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: uent, vent
535    REAL, DIMENSION(nloc, na, na), INTENT(OUT)       :: Sigij, elij
536    REAL, DIMENSION(nloc, na), INTENT(OUT)           :: supmax            ! Highest mixing fraction of mixed
537    ! updraughts with the sign of (h-hp)
538    REAL, DIMENSION(nloc, nd, nd, ntra), INTENT(OUT) :: traent
539    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: Ments, Qents
540    REAL, DIMENSION(nloc, nd, nd), INTENT(OUT)       :: hent
541    INTEGER, DIMENSION(nloc, nd), INTENT(OUT)        :: nent
542
543    !local variables:
544    INTEGER i, j, k, il, im, jm
545    INTEGER num1, num2
546    REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
547    REAL                               :: alt, delp, delm
548    REAL, DIMENSION(nloc)             :: Qmixmax, Rmixmax, sqmrmax
549    REAL, DIMENSION(nloc)             :: Qmixmin, Rmixmin, sqmrmin
550    REAL, DIMENSION(nloc)             :: signhpmh
551    REAL, DIMENSION(nloc)             :: Sx
552    REAL                               :: Scrit2
553    REAL, DIMENSION(nloc)             :: Smid, Sjmin, Sjmax
554    REAL, DIMENSION(nloc)             :: Sbef, sup, smin
555    REAL, DIMENSION(nloc)             :: ASij, ASij_inv, smax, Scrit
556    REAL, DIMENSION(nloc, nd, nd)     :: Sij
557    REAL, DIMENSION(nloc, nd)         :: csum
558    REAL                               :: awat
559    REAL                               :: cpm         !Mixed draught heat capacity
560    REAL                               :: Tm          !Mixed draught temperature
561    LOGICAL, DIMENSION(nloc)          :: lwork
562
563    REAL amxupcrit, df, ff
564    INTEGER nstep
565
566    INTEGER, SAVE                                       :: igout = 1
567    !$omp THREADPRIVATE(igout)
568
569    ! --   Mixing probability distribution functions
570
571    REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
572
573    Qcoef1(F) = tanh(F/gammas)
574    Qcoef2(F) = (tanh(F/gammas) + gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
575    QFF(F) = max(min(F, 1.), 0.)
576    QFFf(F) = min(QFF(F), scut)
577    Qmix1(F) = (tanh((QFF(F) - Fmax)/gammas) + Qcoef1max)/Qcoef2max
578    Rmix1(F) = (gammas*log(cosh((QFF(F) - Fmax)/gammas)) + QFF(F)*Qcoef1max)/Qcoef2max
579    Qmix2(F) = -log(1.-QFFf(F))/scut
580    Rmix2(F) = (QFFf(F) + (1.-QFF(F))*log(1.-QFFf(F)))/scut
581    Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
582    Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
583
584    INTEGER, SAVE :: ifrst
585    DATA ifrst/0/
586    !$omp THREADPRIVATE(ifrst)
587
588    ! =====================================================================
589    ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
590    ! =====================================================================
591
592    ! -- Initialize mixing PDF coefficients
593    IF (ifrst == 0) THEN
594      ifrst = 1
595      Qcoef1max = Qcoef1(Fmax)
596      Qcoef2max = Qcoef2(Fmax)
597      !<jyg
598      print *, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
599        fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
600      !>jyg
601      !
602    END IF
603
604    ! ori        do 360 i=1,ncum*nlp
605    DO j = 1, nl
606    DO i = 1, ncum
607      nent(i, j) = 0
608      ! in convect3, m is computed in cv3_closure
609      ! ori          m(i,1)=0.0
610    END DO
611    END DO
612
613    ! ori      do 400 k=1,nlp
614    ! ori       do 390 j=1,nlp
615    DO j = 1, nl
616    DO k = 1, nl
617    DO i = 1, ncum
618      Qent(i, k, j) = rr(i, j)
619      uent(i, k, j) = u(i, j)
620      vent(i, k, j) = v(i, j)
621      elij(i, k, j) = 0.0
622      hent(i, k, j) = 0.0
623      !AC !            Ment(i,k,j)=0.0
624      !AC !            Sij(i,k,j)=0.0
625    END DO
626    END DO
627    END DO
628
629    !AC !
630    Ment(1:ncum, 1:nd, 1:nd) = 0.0
631    Sij(1:ncum, 1:nd, 1:nd) = 0.0
632    !AC !
633    !ym
634    Sigij(1:ncum, 1:nd, 1:nd) = 0.0
635    !ym
636
637    !jyg !  DO k = 1, ntra
638    !jyg !    DO j = 1, nd  ! instead nlp
639    !jyg !      DO i = 1, nd  ! instead nlp
640    !jyg !        DO il = 1, ncum
641    !jyg !          traent(il, i, j, k) = tra(il, j, k)
642    !jyg !        END DO
643    !jyg !      END DO
644    !jyg !    END DO
645    !jyg !  END DO
646
647    ! =====================================================================
648    ! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
649    ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
650    ! --- FRACTION (Sij)
651    ! =====================================================================
652
653    DO i = minorig + 1, nl
654
655      IF (ok_entrain) THEN
656      DO j = minorig, nl
657      DO il = 1, ncum
658      IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (j >= (icb(il) - 1)) &
659          .AND. (j <= inb(il))) THEN
660
661        ! !            rti = qnk(il) - ep(il, i)*clw(il, i)
662        rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
663        bf2 = 1.+lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd)
664        !jyg(from aj)<
665        IF (cvflag_ice) THEN
666          ! print*,cvflag_ice,'cvflag_ice dans do 700'
667          IF (t(il, j) <= 263.15) THEN
668            bf2 = 1.+(lf(il, j) + lv(il, j))*(lv(il, j) + frac(il, j)* &
669                                              lf(il, j))*rs(il, j)/(rrv*t(il, j)*t(il, j)*cpd)
670          END IF
671        END IF
672        !>jyg
673        anum = h(il, j) - hp(il, i) + (cpv - cpd)*t(il, j)*(rti - rr(il, j))
674        denom = h(il, i) - hp(il, i) + (cpd - cpv)*(rr(il, i) - rti)*t(il, j)
675        dei = denom
676        IF (abs(dei) < 0.01) dei = 0.01
677        Sij(il, i, j) = anum/dei
678        Sij(il, i, i) = 1.0
679        altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j)
680        altem = altem/bf2
681        cwat = clw(il, j)*(1.-ep(il, j))
682        stemp = Sij(il, i, j)
683        IF ((stemp < 0.0 .OR. stemp > 1.0 .OR. altem > cwat) .AND. j > i) THEN
684          !jyg(from aj)<
685          IF (cvflag_ice) THEN
686            anum = anum - (lv(il, j) + frac(il, j)*lf(il, j))*(rti - rs(il, j) - cwat*bf2)
687            denom = denom + (lv(il, j) + frac(il, j)*lf(il, j))*(rr(il, i) - rti)
688          ELSE
689            anum = anum - lv(il, j)*(rti - rs(il, j) - cwat*bf2)
690            denom = denom + lv(il, j)*(rr(il, i) - rti)
691          END IF
692          !>jyg
693          IF (abs(denom) < 0.01) denom = 0.01
694          Sij(il, i, j) = anum/denom
695          altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il, i, j))*rti - rs(il, j)
696          altem = altem - (bf2 - 1.)*cwat
697        END IF
698        IF (Sij(il, i, j) > 0.0) THEN
699          ! ! !                 Ment(il,i,j)=m(il,i)
700          Ment(il, i, j) = 1.
701          elij(il, i, j) = altem
702          elij(il, i, j) = amax1(0.0, elij(il, i, j))
703          nent(il, i) = nent(il, i) + 1
704        END IF
705
706        Sij(il, i, j) = amax1(0.0, Sij(il, i, j))
707        Sij(il, i, j) = amin1(1.0, Sij(il, i, j))
708      ELSE IF (j > i) THEN
709        IF (prt_level >= 10) THEN
710          print *, 'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il, i, j)
711        ENDIF
712      END IF  ! new
713      END DO
714      END DO
715      ELSE   ! (ok_entrain)
716      DO il = 1, ncum
717        nent(il, i) = 0
718      ENDDO
719      ENDIF  ! (ok_entrain)
720
721      !jygdebug<
722      IF (prt_level >= 10) THEN
723        print *, 'cv3p_mixing i, nent(i), icb, inb ', i, nent(igout, i), icb(igout), inb(igout)
724        IF (nent(igout, i) .gt. 0) THEN
725          print *, 'i,(j,Sij(i,j),j=icb-1,inb) ', i, (j, Sij(igout, i, j), j=icb(igout) - 1, inb(igout))
726        ENDIF
727      ENDIF
728      !>jygdebug
729
730      ! ***   if no air can entrain at level i assume that updraft detrains  ***
731      ! ***   at that level and calculate detrained air flux and properties  ***
732
733      ! @      do 170 i=icb(il),inb(il)
734
735      DO il = 1, ncum
736      IF ((i >= icb(il)) .AND. (i <= inb(il)) .AND. (nent(il, i) == 0)) THEN
737        ! @      if(nent(il,i).eq.0)then
738        ! ! !       Ment(il,i,i)=m(il,i)
739        Ment(il, i, i) = 1.
740        ! !        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
741        Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
742        uent(il, i, i) = unk(il)
743        vent(il, i, i) = vnk(il)
744        IF (fl_cor_ebil .GE. 2) THEN
745          hent(il, i, i) = hp(il, i)
746        ENDIF
747        elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
748        Sij(il, i, i) = 0.0
749      END IF
750      END DO
751    END DO  ! i = minorig + 1, nl
752
753    !jyg !  DO j = 1, ntra
754    !jyg !    DO i = minorig + 1, nl
755    !jyg !      DO il = 1, ncum
756    !jyg !        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
757    !jyg !          traent(il, i, i, j) = tra(il, nk(il), j)
758    !jyg !        END IF
759    !jyg !      END DO
760    !jyg !    END DO
761    !jyg !  END DO
762
763    DO j = minorig, nl
764    DO i = minorig, nl
765    DO il = 1, ncum
766    IF ((j >= (icb(il) - 1)) .AND. (j <= inb(il)) .AND. &
767        (i >= icb(il)) .AND. (i <= inb(il))) THEN
768      Sigij(il, i, j) = Sij(il, i, j)
769    END IF
770    END DO
771    END DO
772    END DO
773    ! @      enddo
774
775    ! @170   continue
776
777    ! =====================================================================
778    ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
779    ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
780    ! =====================================================================
781
782    CALL zilch(csum, nloc*nd)
783
784    DO il = 1, ncum
785      lwork(il) = .FALSE.
786    END DO
787
788    ! ---------------------------------------------------------------
789    DO i = minorig + 1, nl       !Loop on origin level "i"
790      ! ---------------------------------------------------------------
791
792      num1 = 0
793      DO il = 1, ncum
794        IF (i >= icb(il) .AND. i <= inb(il)) num1 = num1 + 1
795      END DO
796      IF (num1 <= 0) GO TO 789
797
798      !JYG1    Find maximum of SIJ for J>I, if any.
799
800      Sx(:) = 0.
801
802      DO il = 1, ncum
803      IF (i >= icb(il) .AND. i <= inb(il)) THEN
804        signhpmh(il) = sign(1., hp(il, i) - h(il, i))
805        Sbef(il) = max(0., signhpmh(il))
806      END IF
807      END DO
808
809      DO j = i + 1, nl
810      DO il = 1, ncum
811      IF (i >= icb(il) .AND. i <= inb(il) .AND. j <= inb(il)) THEN
812      IF (Sbef(il) < Sij(il, i, j)) THEN
813        Sx(il) = max(Sij(il, i, j), Sx(il))
814      END IF
815      Sbef(il) = Sij(il, i, j)
816      END IF
817      END DO
818      END DO
819
820      DO il = 1, ncum
821      IF (i >= icb(il) .AND. i <= inb(il)) THEN
822        lwork(il) = (nent(il, i) /= 0)
823        ! !        rti = qnk(il) - ep(il, i)*clw(il, i)
824        rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
825        !jyg<
826        IF (cvflag_ice) THEN
827
828          anum = h(il, i) - hp(il, i) - (lv(il, i) + frac(il, i)*lf(il, i))* &
829                 (rti - rs(il, i)) + (cpv - cpd)*t(il, i)*(rti - rr(il, i))
830          denom = h(il, i) - hp(il, i) + (lv(il, i) + frac(il, i)*lf(il, i))* &
831                  (rr(il, i) - rti) + (cpd - cpv)*t(il, i)*(rr(il, i) - rti)
832        ELSE
833
834          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti - rs(il, i)) + &
835                 (cpv - cpd)*t(il, i)*(rti - rr(il, i))
836          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il, i) - rti) + &
837                  (cpd - cpv)*t(il, i)*(rr(il, i) - rti)
838        END IF
839        !>jyg
840        IF (abs(denom) < 0.01) denom = 0.01
841        Scrit(il) = min(anum/denom, 1.)
842        alt = rti - rs(il, i) + Scrit(il)*(rr(il, i) - rti)
843
844        !JYG1    Find new critical value Scrit2
845        !         such that : Sij > Scrit2  => mixed draught will detrain at J<I
846        !                     Sij < Scrit2  => mixed draught will detrain at J>I
847
848        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
849                 Scrit(il)*max(0., signhpmh(il))
850
851        Scrit(il) = Scrit2
852
853        !JYG    Correction pour la nouvelle logique; la correction pour ALT
854        ! est un peu au hazard
855        IF (Scrit(il) <= 0.0) Scrit(il) = 0.0
856        IF (alt <= 0.0) Scrit(il) = 1.0
857
858        smax(il) = 0.0
859        ASij(il) = 0.0
860        sup(il) = 0.       ! upper S-value reached by descending draughts
861      END IF
862      END DO
863
864      ! ---------------------------------------------------------------
865      DO j = minorig, nl          !Loop on destination level "j"
866        ! ---------------------------------------------------------------
867
868        num2 = 0
869        DO il = 1, ncum
870          IF (i >= icb(il) .AND. i <= inb(il) .AND. &
871              j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
872              lwork(il)) num2 = num2 + 1
873        END DO
874        IF (num2 <= 0) GO TO 175
875
876        ! -----------------------------------------------
877        IF (j > i) THEN
878          ! -----------------------------------------------
879          DO il = 1, ncum
880          IF (i >= icb(il) .AND. i <= inb(il) .AND. &
881              j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
882              lwork(il)) THEN
883          IF (Sij(il, i, j) > 0.0) THEN
884            Smid(il) = min(Sij(il, i, j), Scrit(il))
885            Sjmax(il) = Smid(il)
886            Sjmin(il) = Smid(il)
887            IF (Smid(il) < smin(il) .AND. Sij(il, i, j + 1) < Smid(il)) THEN
888              smin(il) = Smid(il)
889              Sjmax(il) = min((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j), Scrit(il))
890              Sjmin(il) = max((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
891              Sjmin(il) = min(Sjmin(il), Scrit(il))
892              Sbef(il) = Sij(il, i, j)
893            END IF
894          END IF
895          END IF
896          END DO
897          ! -----------------------------------------------
898        ELSE IF (j == i) THEN
899          ! -----------------------------------------------
900          DO il = 1, ncum
901          IF (i >= icb(il) .AND. i <= inb(il) .AND. &
902              j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
903              lwork(il)) THEN
904          IF (Sij(il, i, j) > 0.0) THEN
905            Smid(il) = 1.
906            Sjmin(il) = max((Sij(il, i, j - 1) + Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
907                        min((Sij(il, i, j + 1) + Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
908            Sjmin(il) = max(Sjmin(il), sup(il))
909            Sjmax(il) = 1.
910
911            ! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
912            Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
913
914            smin(il) = 1.
915            Sbef(il) = max(0., signhpmh(il))
916            supmax(il, i) = sign(Scrit(il), -signhpmh(il))
917          END IF
918          END IF
919          END DO
920          ! -----------------------------------------------
921        ELSE IF (j < i) THEN
922          ! -----------------------------------------------
923          DO il = 1, ncum
924          IF (i >= icb(il) .AND. i <= inb(il) .AND. &
925              j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
926              lwork(il)) THEN
927          IF (Sij(il, i, j) > 0.0) THEN
928            Smid(il) = max(Sij(il, i, j), Scrit(il))
929            Sjmax(il) = Smid(il)
930            Sjmin(il) = Smid(il)
931            IF (Smid(il) > smax(il) .AND. Sij(il, i, j + 1) > Smid(il)) THEN
932              smax(il) = Smid(il)
933              Sjmax(il) = max((Sij(il, i, j + 1) + Sij(il, i, j))/2., Sij(il, i, j))
934              Sjmax(il) = max(Sjmax(il), Scrit(il))
935              Sjmin(il) = min((Sbef(il) + Sij(il, i, j))/2., Sij(il, i, j))
936              Sjmin(il) = max(Sjmin(il), Scrit(il))
937              Sbef(il) = Sij(il, i, j)
938            END IF
939            IF (abs(Sjmin(il) - Sjmax(il)) > 1.E-10) &
940              sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
941          END IF
942          END IF
943          END DO
944          ! -----------------------------------------------
945        END IF
946        ! -----------------------------------------------
947
948        DO il = 1, ncum
949        IF (i >= icb(il) .AND. i <= inb(il) .AND. &
950            j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
951            lwork(il)) THEN
952        IF (Sij(il, i, j) > 0.0) THEN
953          ! !            rti = qnk(il) - ep(il, i)*clw(il, i)
954          rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
955          Qmixmax(il) = Qmix(Sjmax(il))
956          Qmixmin(il) = Qmix(Sjmin(il))
957          Rmixmax(il) = Rmix(Sjmax(il))
958          Rmixmin(il) = Rmix(Sjmin(il))
959          sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
960          sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
961
962          Ment(il, i, j) = abs(Qmixmax(il) - Qmixmin(il))*Ment(il, i, j)
963
964          ! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
965          IF (abs(Qmixmax(il) - Qmixmin(il)) > 1.E-10) THEN
966            Sigij(il, i, j) = (sqmrmax(il) - sqmrmin(il))/(Qmixmax(il) - Qmixmin(il))
967          ELSE
968            Sigij(il, i, j) = 0.
969          END IF
970
971          ! --    Compute Qent, uent, vent according to the true mixing fraction
972          Qent(il, i, j) = (1.-Sigij(il, i, j))*rti + Sigij(il, i, j)*rr(il, i)
973          uent(il, i, j) = (1.-Sigij(il, i, j))*unk(il) + Sigij(il, i, j)*u(il, i)
974          vent(il, i, j) = (1.-Sigij(il, i, j))*vnk(il) + Sigij(il, i, j)*v(il, i)
975
976          ! --     Compute liquid water static energy of mixed draughts
977          !    IF (j .GT. i) THEN
978          !      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
979          !      awat=amax1(awat,0.0)
980          !    ELSE
981          !      awat = 0.
982          !    ENDIF
983          !    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
984          !    :         + Sigij(il,i,j)*H(il,i)
985          !    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
986          !IM 301008 beg
987          hent(il, i, j) = (1.-Sigij(il, i, j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
988
989          !jyg<
990          !            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
991          !            elij(il, i, j) = elij(il, i, j) + &
992          !                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
993          !                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
994          !            elij(il, i, j) = elij(il, i, j) / &
995          !                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
996          !                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
997          !
998          !       Computation of condensate amount Elij, taking into account the ice fraction frac
999          !       Warning : the same saturation humidity rs is used over both liquid water and ice; this
1000          !                 should be corrected.
1001          !
1002          !  Heat capacity of mixed draught
1003          cpm = cpd + Qent(il, i, j)*(cpv - cpd)
1004          !
1005          IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
1006            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
1007            elij(il, i, j) = elij(il, i, j) + &
1008                             (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
1009                             rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
1010            elij(il, i, j) = elij(il, i, j)/ &
1011                             (1.+(lv(il, j) + frac(il, j)*lf(il, j))*lv(il, j)*rs(il, j)/ &
1012                              (cpm*rrv*t(il, j)*t(il, j)))
1013          ELSE
1014            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
1015            elij(il, i, j) = elij(il, i, j) + &
1016                             (h(il, j) - hent(il, i, j) + (cpv - cpd)*(Qent(il, i, j) - rr(il, j))*t(il, j))* &
1017                             rs(il, j)*lv(il, j)/(cpm*rrv*t(il, j)*t(il, j))
1018            elij(il, i, j) = elij(il, i, j)/ &
1019                             (1.+lv(il, j)*lv(il, j)*rs(il, j)/ &
1020                              (cpm*rrv*t(il, j)*t(il, j)))
1021          ENDIF
1022          !>jyg
1023          elij(il, i, j) = max(elij(il, i, j), 0.)
1024
1025          elij(il, i, j) = min(elij(il, i, j), Qent(il, i, j))
1026
1027          IF (j > i) THEN
1028            awat = elij(il, i, j) - (1.-ep(il, j))*clw(il, j)
1029            awat = amax1(awat, 0.0)
1030          ELSE
1031            awat = 0.
1032          END IF
1033
1034          ! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
1035          ! :         t(il,j))
1036
1037          !jyg<
1038          !            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
1039          ! Mixed draught temperature at level j
1040          IF (cvflag_ice .and. frac(il, j) .gt. 0.) THEN
1041            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))
1042            hent(il, i, j) = hent(il, i, j) + (lv(il, j) + frac(il, j)*lf(il, j) + (cpd - cpv)*Tm)*awat
1043          ELSE
1044            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))
1045            hent(il, i, j) = hent(il, i, j) + (lv(il, j) + (cpd - cpv)*Tm)*awat
1046          ENDIF
1047          !>jyg
1048
1049          !IM 301008 end
1050
1051          ! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
1052          ! :               i,j,hent(il,i,j),Sigij(il,i,j)
1053
1054          ! --      ASij is the integral of P(F) over the relevant F interval
1055          ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
1056                                    Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
1057
1058        END IF
1059        END IF
1060        END DO
1061        !jyg !      DO k = 1, ntra
1062        !jyg !        DO il = 1, ncum
1063        !jyg !          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1064        !jyg !              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1065        !jyg !              lwork(il)) THEN
1066        !jyg !            IF (Sij(il,i,j)>0.0) THEN
1067        !jyg !              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
1068        !jyg !                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
1069        !jyg !            END IF
1070        !jyg !          END IF
1071        !jyg !        END DO
1072        !jyg !      END DO
1073
1074        ! --    If I=J (detrainement and entrainement at the same level), then only the
1075        ! --    adiabatic ascent part of the mixture is considered
1076        IF (i == j) THEN
1077        DO il = 1, ncum
1078        IF (i >= icb(il) .AND. i <= inb(il) .AND. &
1079            j >= (icb(il) - 1) .AND. j <= inb(il) .AND. &
1080            lwork(il)) THEN
1081        IF (Sij(il, i, j) > 0.0) THEN
1082          ! !              rti = qnk(il) - ep(il, i)*clw(il, i)
1083          rti = qta(il, i - 1) - ep(il, i)*clw(il, i)
1084          ! ! !             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
1085          Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il)) + Rmixmax(il) - &
1086                               Qmixmin(il)*(1.-Sjmin(il)) - Rmixmin(il))
1087          Qent(il, i, i) = rti
1088          uent(il, i, i) = unk(il)
1089          vent(il, i, i) = vnk(il)
1090          hent(il, i, i) = hp(il, i)
1091          elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
1092          Sigij(il, i, i) = 0.
1093        END IF
1094        END IF
1095        END DO
1096        !jyg !        DO k = 1, ntra
1097        !jyg !          DO il = 1, ncum
1098        !jyg !            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1099        !jyg !                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1100        !jyg !                lwork(il)) THEN
1101        !jyg !              IF (Sij(il,i,j)>0.0) THEN
1102        !jyg !                traent(il, i, i, k) = tra(il, nk(il), k)
1103        !jyg !              END IF
1104        !jyg !            END IF
1105        !jyg !          END DO
1106        !jyg !        END DO
1107
1108        END IF
1109
1110        ! ---------------------------------------------------------------
1111175   END DO         ! End loop on destination level "j"
1112      ! ---------------------------------------------------------------
1113
1114      DO il = 1, ncum
1115      IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il)) THEN
1116        ASij(il) = amax1(1.0E-16, ASij(il))
1117        !jyg+lluis<
1118        ! !        ASij(il) = 1.0/ASij(il)
1119        ASij_inv(il) = 1.0/ASij(il)
1120        !   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
1121        IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
1122        !>jyg+lluis
1123        csum(il, i) = 0.0
1124      END IF
1125      END DO
1126
1127      DO j = minorig, nl
1128      DO il = 1, ncum
1129      IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
1130          j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
1131        !jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
1132        Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
1133      END IF
1134      END DO
1135      END DO
1136
1137      DO j = minorig, nl
1138      DO il = 1, ncum
1139      IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. &
1140          j >= (icb(il) - 1) .AND. j <= inb(il)) THEN
1141        csum(il, i) = csum(il, i) + Ment(il, i, j)
1142      END IF
1143      END DO
1144      END DO
1145
1146      DO il = 1, ncum
1147      IF (i >= icb(il) .AND. i <= inb(il) .AND. lwork(il) .AND. csum(il, i) < 1.) THEN
1148        ! cc     :     .and. csum(il,i).lt.m(il,i) ) then
1149        nent(il, i) = 0
1150        ! cc        Ment(il,i,i)=m(il,i)
1151        Ment(il, i, i) = 1.
1152        ! !        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1153        Qent(il, i, i) = qta(il, i - 1) - ep(il, i)*clw(il, i)
1154        uent(il, i, i) = unk(il)
1155        vent(il, i, i) = vnk(il)
1156        elij(il, i, i) = clw(il, i)*(1.-ep(il, i))
1157        IF (fl_cor_ebil .GE. 2) THEN
1158          hent(il, i, i) = hp(il, i)
1159          Sigij(il, i, i) = 0.0
1160        ELSE
1161          Sij(il, i, i) = 0.0
1162        ENDIF
1163      END IF
1164      END DO  ! il
1165
1166      !jyg !    DO j = 1, ntra
1167      !jyg !      DO il = 1, ncum
1168      !jyg !        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1169      !jyg ! ! cc     :     .and. csum(il,i).lt.m(il,i) ) then
1170      !jyg !          traent(il, i, i, j) = tra(il, nk(il), j)
1171      !jyg !        END IF
1172      !jyg !      END DO
1173      !jyg !    END DO
1174
1175      ! ---------------------------------------------------------------
1176789 END DO               ! End loop on origin level "i"
1177    ! ---------------------------------------------------------------
1178
1179    RETURN
1180  END SUBROUTINE
1181end module
Note: See TracBrowser for help on using the repository browser.