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

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

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