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

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

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