source: LMDZ6/branches/cirrus/libf/phylmd/cv3p_mixing.F90 @ 5185

Last change on this file since 5185 was 3496, checked in by jyg, 6 years ago

Implementation of the ejection of liquid precipitation from the adiabatic ascents.
New flags:
+cvflag_prec_eject: logical

n -> old code, y -> new code

+ejectliq: real; possible values 0. & 1.

  1. -> no liquid precipitation is ejected
  2. -> all liquid precipitation is ejected

+ejectice: real; any value between 0. and 1.

fraction of solid precipitation ejected at each level

Note that the adiabatic ascent mass flux decrease due to precipitation ejection is not taken into account.

Attempts to do it led to water conservation violation.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.3 KB
Line 
1SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
3                       unk, vnk, hp, tv, tvp, ep, clw, sig, &
4                       Ment, Qent, hent, uent, vent, nent, &
5                       Sigij, elij, supmax, Ments, Qents, traent)
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
15  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
16  USE ioipsl_getin_p_mod, ONLY: getin_p
17  USE add_phys_tend_mod, ONLY: fl_cor_ebil
18
19  IMPLICIT NONE
20
21  include "cvthermo.h"
22  include "cv3param.h"
23  include "YOMCST2.h"
24  include "cvflag.h"
25
26!inputs:
27  INTEGER, INTENT (IN)                               :: ncum, nd, na
28  INTEGER, INTENT (IN)                               :: ntra, nloc
29  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
30  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
31  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
32  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
33  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
34  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
35  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
36  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
37  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
38  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
39  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
40  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
41  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
42  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
43  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
44
45!outputs:
46  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
47  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
48  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
49  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax           ! Highest mixing fraction of mixed
50                                                                         ! updraughts with the sign of (h-hp)
51  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
52  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
53  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
54  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
55
56!local variables:
57  INTEGER i, j, k, il, im, jm
58  INTEGER num1, num2
59  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
60  REAL                               :: alt, delp, delm
61  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
62  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
63  REAL, DIMENSION (nloc)             :: signhpmh
64  REAL, DIMENSION (nloc)             :: Sx
65  REAL                               :: Scrit2
66  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
67  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
68  REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
69  REAL, DIMENSION (nloc, nd, nd)     :: Sij
70  REAL, DIMENSION (nloc, nd)         :: csum
71  REAL                               :: awat
72  REAL                               :: cpm        !Mixed draught heat capacity
73  REAL                               :: Tm         !Mixed draught temperature
74  LOGICAL, DIMENSION (nloc)          :: lwork
75
76  REAL amxupcrit, df, ff
77  INTEGER nstep
78
79  INTEGER,SAVE                                       :: igout=1
80!$OMP THREADPRIVATE(igout)
81
82! --   Mixing probability distribution functions
83
84  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
85
86  Qcoef1(F) = tanh(F/gammas)
87  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
88  QFF(F) = max(min(F,1.), 0.)
89  QFFf(F) = min(QFF(F), scut)
90  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
91  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
92  Qmix2(F) = -log(1.-QFFf(F))/scut
93  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
94  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
95  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
96
97  INTEGER, SAVE :: ifrst
98  DATA ifrst/0/
99!$OMP THREADPRIVATE(ifrst)
100
101
102! =====================================================================
103! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
104! =====================================================================
105
106! -- Initialize mixing PDF coefficients
107  IF (ifrst==0) THEN
108    ifrst = 1
109    Qcoef1max = Qcoef1(Fmax)
110    Qcoef2max = Qcoef2(Fmax)
111!<jyg
112   print*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
113            fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
114!>jyg
115!
116  END IF
117
118
119! ori        do 360 i=1,ncum*nlp
120  DO j = 1, nl
121    DO i = 1, ncum
122      nent(i, j) = 0
123! in convect3, m is computed in cv3_closure
124! ori          m(i,1)=0.0
125    END DO
126  END DO
127
128! ori      do 400 k=1,nlp
129! ori       do 390 j=1,nlp
130  DO j = 1, nl
131    DO k = 1, nl
132      DO i = 1, ncum
133        Qent(i, k, j) = rr(i, j)
134        uent(i, k, j) = u(i, j)
135        vent(i, k, j) = v(i, j)
136        elij(i, k, j) = 0.0
137        hent(i, k, j) = 0.0
138!AC!            Ment(i,k,j)=0.0
139!AC!            Sij(i,k,j)=0.0
140      END DO
141    END DO
142  END DO
143
144!AC!
145  Ment(1:ncum, 1:nd, 1:nd) = 0.0
146  Sij(1:ncum, 1:nd, 1:nd) = 0.0
147!AC!
148!ym
149  Sigij(1:ncum, 1:nd, 1:nd) = 0.0
150!ym
151
152!jyg!  DO k = 1, ntra
153!jyg!    DO j = 1, nd ! instead nlp
154!jyg!      DO i = 1, nd ! instead nlp
155!jyg!        DO il = 1, ncum
156!jyg!          traent(il, i, j, k) = tra(il, j, k)
157!jyg!        END DO
158!jyg!      END DO
159!jyg!    END DO
160!jyg!  END DO
161
162! =====================================================================
163! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
164! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
165! --- FRACTION (Sij)
166! =====================================================================
167
168  DO i = minorig + 1, nl
169
170    IF (ok_entrain) THEN
171      DO j = minorig, nl
172        DO il = 1, ncum
173          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
174                           .AND. (j<=inb(il))) THEN
175
176!!            rti = qnk(il) - ep(il, i)*clw(il, i)
177            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
178            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
179!jyg(from aj)<
180            IF (cvflag_ice) THEN
181! print*,cvflag_ice,'cvflag_ice dans do 700'
182              IF (t(il,j)<=263.15) THEN
183                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
184                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
185              END IF
186            END IF
187!>jyg
188            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
189            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
190            dei = denom
191            IF (abs(dei)<0.01) dei = 0.01
192            Sij(il, i, j) = anum/dei
193            Sij(il, i, i) = 1.0
194            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
195            altem = altem/bf2
196            cwat = clw(il, j)*(1.-ep(il,j))
197            stemp = Sij(il, i, j)
198            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
199!jyg(from aj)<
200              IF (cvflag_ice) THEN
201                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
202                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
203              ELSE
204                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
205                denom = denom + lv(il, j)*(rr(il,i)-rti)
206              END IF
207!>jyg
208              IF (abs(denom)<0.01) denom = 0.01
209              Sij(il, i, j) = anum/denom
210              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
211              altem = altem - (bf2-1.)*cwat
212            END IF
213            IF (Sij(il,i,j)>0.0) THEN
214!!!                 Ment(il,i,j)=m(il,i)
215              Ment(il, i, j) = 1.
216              elij(il, i, j) = altem
217              elij(il, i, j) = amax1(0.0, elij(il,i,j))
218              nent(il, i) = nent(il, i) + 1
219            END IF
220
221            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
222            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
223          ELSE IF (j > i) THEN
224            IF (prt_level >= 10) THEN
225              print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
226            ENDIF
227          END IF ! new
228        END DO
229      END DO
230    ELSE  ! (ok_entrain)
231      DO il = 1,ncum
232        nent(il,i) = 0
233      ENDDO
234    ENDIF ! (ok_entrain)
235
236!jygdebug<
237    IF (prt_level >= 10) THEN
238      print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
239      IF (nent(igout,i) .gt. 0) THEN
240        print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
241      ENDIF
242    ENDIF
243!>jygdebug
244
245! ***   if no air can entrain at level i assume that updraft detrains  ***
246! ***   at that level and calculate detrained air flux and properties  ***
247
248
249! @      do 170 i=icb(il),inb(il)
250
251    DO il = 1, ncum
252      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
253! @      if(nent(il,i).eq.0)then
254!!!       Ment(il,i,i)=m(il,i)
255        Ment(il, i, i) = 1.
256!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
257        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
258        uent(il, i, i) = unk(il)
259        vent(il, i, i) = vnk(il)
260        IF (fl_cor_ebil .GE. 2) THEN
261          hent(il, i, i) = hp(il,i)
262        ENDIF
263        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
264        Sij(il, i, i) = 0.0
265      END IF
266    END DO
267  END DO ! i = minorig + 1, nl
268
269!jyg!  DO j = 1, ntra
270!jyg!    DO i = minorig + 1, nl
271!jyg!      DO il = 1, ncum
272!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
273!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
274!jyg!        END IF
275!jyg!      END DO
276!jyg!    END DO
277!jyg!  END DO
278
279  DO j = minorig, nl
280    DO i = minorig, nl
281      DO il = 1, ncum
282        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
283            (i>=icb(il)) .AND. (i<=inb(il))) THEN
284          Sigij(il, i, j) = Sij(il, i, j)
285        END IF
286      END DO
287    END DO
288  END DO
289! @      enddo
290
291! @170   continue
292
293! =====================================================================
294! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
295! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
296! =====================================================================
297
298  CALL zilch(csum, nloc*nd)
299
300  DO il = 1, ncum
301    lwork(il) = .FALSE.
302  END DO
303
304! ---------------------------------------------------------------
305  DO i = minorig + 1, nl      !Loop on origin level "i"
306! ---------------------------------------------------------------
307
308    num1 = 0
309    DO il = 1, ncum
310      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
311    END DO
312    IF (num1<=0) GO TO 789
313
314
315!JYG1    Find maximum of SIJ for J>I, if any.
316
317    Sx(:) = 0.
318
319    DO il = 1, ncum
320      IF (i>=icb(il) .AND. i<=inb(il)) THEN
321        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
322        Sbef(il) = max(0., signhpmh(il))
323      END IF
324    END DO
325
326    DO j = i + 1, nl
327      DO il = 1, ncum
328        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
329          IF (Sbef(il)<Sij(il,i,j)) THEN
330            Sx(il) = max(Sij(il,i,j), Sx(il))
331          END IF
332          Sbef(il) = Sij(il, i, j)
333        END IF
334      END DO
335    END DO
336
337
338    DO il = 1, ncum
339      IF (i>=icb(il) .AND. i<=inb(il)) THEN
340        lwork(il) = (nent(il,i)/=0)
341!!        rti = qnk(il) - ep(il, i)*clw(il, i)
342        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
343!jyg<
344        IF (cvflag_ice) THEN
345
346          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
347                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
348          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
349                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
350        ELSE
351
352          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
353                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
354          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
355                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
356        END IF
357!>jyg
358        IF (abs(denom)<0.01) denom = 0.01
359        Scrit(il) = min(anum/denom, 1.)
360        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
361
362!JYG1    Find new critical value Scrit2
363!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
364!                     Sij < Scrit2  => mixed draught will detrain at J>I
365
366        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
367                 Scrit(il)*max(0., signhpmh(il))
368
369        Scrit(il) = Scrit2
370
371!JYG    Correction pour la nouvelle logique; la correction pour ALT
372! est un peu au hazard
373        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
374        IF (alt<=0.0) Scrit(il) = 1.0
375
376        smax(il) = 0.0
377        ASij(il) = 0.0
378        sup(il) = 0.      ! upper S-value reached by descending draughts
379      END IF
380    END DO
381
382! ---------------------------------------------------------------
383    DO j = minorig, nl         !Loop on destination level "j"
384! ---------------------------------------------------------------
385
386      num2 = 0
387      DO il = 1, ncum
388        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
389            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
390            lwork(il)) num2 = num2 + 1
391      END DO
392      IF (num2<=0) GO TO 175
393
394! -----------------------------------------------
395      IF (j>i) THEN
396! -----------------------------------------------
397        DO il = 1, ncum
398          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
399              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
400              lwork(il)) THEN
401            IF (Sij(il,i,j)>0.0) THEN
402              Smid(il) = min(Sij(il,i,j), Scrit(il))
403              Sjmax(il) = Smid(il)
404              Sjmin(il) = Smid(il)
405              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
406                smin(il) = Smid(il)
407                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
408                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
409                Sjmin(il) = min(Sjmin(il), Scrit(il))
410                Sbef(il) = Sij(il, i, j)
411              END IF
412            END IF
413          END IF
414        END DO
415! -----------------------------------------------
416      ELSE IF (j==i) THEN
417! -----------------------------------------------
418        DO il = 1, ncum
419          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
420              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
421              lwork(il)) THEN
422            IF (Sij(il,i,j)>0.0) THEN
423              Smid(il) = 1.
424              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
425                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
426              Sjmin(il) = max(Sjmin(il), sup(il))
427              Sjmax(il) = 1.
428
429! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
430              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
431
432              smin(il) = 1.
433              Sbef(il) = max(0., signhpmh(il))
434              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
435            END IF
436          END IF
437        END DO
438! -----------------------------------------------
439      ELSE IF (j<i) THEN
440! -----------------------------------------------
441        DO il = 1, ncum
442          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
443              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
444              lwork(il)) THEN
445            IF (Sij(il,i,j)>0.0) THEN
446              Smid(il) = max(Sij(il,i,j), Scrit(il))
447              Sjmax(il) = Smid(il)
448              Sjmin(il) = Smid(il)
449              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
450                smax(il) = Smid(il)
451                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
452                Sjmax(il) = max(Sjmax(il), Scrit(il))
453                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
454                Sjmin(il) = max(Sjmin(il), Scrit(il))
455                Sbef(il) = Sij(il, i, j)
456              END IF
457              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
458                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
459            END IF
460          END IF
461        END DO
462! -----------------------------------------------
463      END IF
464! -----------------------------------------------
465
466
467      DO il = 1, ncum
468        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
469            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
470            lwork(il)) THEN
471          IF (Sij(il,i,j)>0.0) THEN
472!!            rti = qnk(il) - ep(il, i)*clw(il, i)
473            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
474            Qmixmax(il) = Qmix(Sjmax(il))
475            Qmixmin(il) = Qmix(Sjmin(il))
476            Rmixmax(il) = Rmix(Sjmax(il))
477            Rmixmin(il) = Rmix(Sjmin(il))
478            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
479            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
480
481            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
482
483! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
484            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
485              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
486            ELSE
487              Sigij(il, i, j) = 0.
488            END IF
489
490! --    Compute Qent, uent, vent according to the true mixing fraction
491            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
492            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
493            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
494
495! --     Compute liquid water static energy of mixed draughts
496!    IF (j .GT. i) THEN
497!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
498!      awat=amax1(awat,0.0)
499!    ELSE
500!      awat = 0.
501!    ENDIF
502!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
503!    :         + Sigij(il,i,j)*H(il,i)
504!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
505!IM 301008 beg
506            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
507
508!jyg<
509!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
510!            elij(il, i, j) = elij(il, i, j) + &
511!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
512!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
513!            elij(il, i, j) = elij(il, i, j) / &
514!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
515!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
516!
517!       Computation of condensate amount Elij, taking into account the ice fraction frac
518!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
519!                 should be corrected.
520!
521!  Heat capacity of mixed draught
522    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
523!
524    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
525            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
526            elij(il, i, j) = elij(il, i, j) + &
527                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
528                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
529            elij(il, i, j) = elij(il, i, j) / &
530                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
531                              (cpm*rrv*t(il,j)*t(il,j)))
532    ELSE
533            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
534            elij(il, i, j) = elij(il, i, j) + &
535                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
536                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
537            elij(il, i, j) = elij(il, i, j) / &
538                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
539                              (cpm*rrv*t(il,j)*t(il,j)))
540    ENDIF
541!>jyg
542            elij(il, i, j) = max(elij(il,i,j), 0.)
543
544            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
545
546            IF (j>i) THEN
547              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
548              awat = amax1(awat, 0.0)
549            ELSE
550              awat = 0.
551            END IF
552
553! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
554! :         t(il,j))
555
556!jyg<
557!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
558! Mixed draught temperature at level j
559    IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
560          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))
561          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
562    ELSE
563          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))
564          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
565    ENDIF
566!>jyg
567
568!IM 301008 end
569
570! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
571! :               i,j,hent(il,i,j),Sigij(il,i,j)
572
573! --      ASij is the integral of P(F) over the relevant F interval
574            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
575                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
576
577          END IF
578        END IF
579      END DO
580!jyg!      DO k = 1, ntra
581!jyg!        DO il = 1, ncum
582!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
583!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
584!jyg!              lwork(il)) THEN
585!jyg!            IF (Sij(il,i,j)>0.0) THEN
586!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
587!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
588!jyg!            END IF
589!jyg!          END IF
590!jyg!        END DO
591!jyg!      END DO
592
593! --    If I=J (detrainement and entrainement at the same level), then only the
594! --    adiabatic ascent part of the mixture is considered
595      IF (i==j) THEN
596        DO il = 1, ncum
597          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
598              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
599              lwork(il)) THEN
600            IF (Sij(il,i,j)>0.0) THEN
601!!              rti = qnk(il) - ep(il, i)*clw(il, i)
602              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
603!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
604              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
605                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
606              Qent(il, i, i) = rti
607              uent(il, i, i) = unk(il)
608              vent(il, i, i) = vnk(il)
609              hent(il, i, i) = hp(il, i)
610              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
611              Sigij(il, i, i) = 0.
612            END IF
613          END IF
614        END DO
615!jyg!        DO k = 1, ntra
616!jyg!          DO il = 1, ncum
617!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
618!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
619!jyg!                lwork(il)) THEN
620!jyg!              IF (Sij(il,i,j)>0.0) THEN
621!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
622!jyg!              END IF
623!jyg!            END IF
624!jyg!          END DO
625!jyg!        END DO
626
627      END IF
628
629! ---------------------------------------------------------------
630175 END DO        ! End loop on destination level "j"
631! ---------------------------------------------------------------
632
633    DO il = 1, ncum
634      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
635        ASij(il) = amax1(1.0E-16, ASij(il))
636!jyg+lluis<
637!!        ASij(il) = 1.0/ASij(il)
638        ASij_inv(il) = 1.0/ASij(il)
639!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
640        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
641!>jyg+lluis
642        csum(il, i) = 0.0
643      END IF
644    END DO
645
646    DO j = minorig, nl
647      DO il = 1, ncum
648        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
649            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
650!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
651          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
652        END IF
653      END DO
654    END DO
655
656    DO j = minorig, nl
657      DO il = 1, ncum
658        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
659            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
660          csum(il, i) = csum(il, i) + Ment(il, i, j)
661        END IF
662      END DO
663    END DO
664
665    DO il = 1, ncum
666      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
667! cc     :     .and. csum(il,i).lt.m(il,i) ) then
668        nent(il, i) = 0
669! cc        Ment(il,i,i)=m(il,i)
670        Ment(il, i, i) = 1.
671!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
672        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
673        uent(il, i, i) = unk(il)
674        vent(il, i, i) = vnk(il)
675        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
676        IF (fl_cor_ebil .GE. 2) THEN
677          hent(il, i, i) = hp(il,i)
678          Sigij(il, i, i) = 0.0
679        ELSE
680          Sij(il, i, i) = 0.0
681        ENDIF
682      END IF
683    END DO ! il
684
685!jyg!    DO j = 1, ntra
686!jyg!      DO il = 1, ncum
687!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
688!jyg!! cc     :     .and. csum(il,i).lt.m(il,i) ) then
689!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
690!jyg!        END IF
691!jyg!      END DO
692!jyg!    END DO
693
694! ---------------------------------------------------------------
695789 END DO              ! End loop on origin level "i"
696! ---------------------------------------------------------------
697
698
699  RETURN
700END SUBROUTINE cv3p_mixing
701
Note: See TracBrowser for help on using the repository browser.