source: LMDZ5/trunk/libf/phylmd/cv3p_mixing.F90 @ 3712

Last change on this file since 3712 was 2905, checked in by musat, 7 years ago

Bug correction on ok_entrain making an MPI error

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