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

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

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