source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/cv3p_mixing.F90 @ 3433

Last change on this file since 3433 was 2478, checked in by fhourdin, 8 years ago

Correction d'un bug de non initialisation d'un tableau qui faisait
planter le modèle en haute résolution sur 8 thread. Bug identifié
par Yann Meurdesoif.

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