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

Last change on this file since 3715 was 3715, checked in by adurocher, 4 years ago

cv3p_mixing_new : change array temporaries into scalars

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