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

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

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