source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/cv3p_mixing.F90 @ 5141

Last change on this file since 5141 was 5141, checked in by abarral, 8 weeks ago

Put cvthermo.h, cv30param.h, cv3param.h into modules

  • Property svn:keywords set to Id
File size: 50.7 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#ifdef ISO
7                         ,xt,xtta,xtclw,xtent,xtelij  &
8#endif         
9                         )
10! **************************************************************
11! *
12! CV3P_MIXING : compute mixed draught properties and,         *
13! within a scaling factor, mixed draught        *
14! mass fluxes.                                  *
15! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
16! modified by :                                               *
17! **************************************************************
18
19  USE lmdz_print_control, ONLY: mydebug=>debug , lunout, prt_level
20  USE lmdz_ioipsl_getin_p, ONLY: getin_p
21  USE add_phys_tend_mod, ONLY: fl_cor_ebil
22#ifdef ISO
23  USE infotrac_phy, ONLY: ntraciso=>ntiso
24  USE isotopes_mod, ONLY: pxtmelt,pxtice
25  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall
26#ifdef ISOTRAC
27  USE infotrac_phy, ONLY: niso
28  USE isotrac_mod, ONLY: option_cond,izone_cond, &
29        option_tmin,option_traceurs,seuil_tag_tmin, &
30        index_zone,index_iso
31  USE isotrac_routines_mod, ONLY: iso_recolorise_condensation
32  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac
33#endif
34#ifdef ISOVERIF
35       USE isotopes_verif_mod
36!, ONLY: iso_verif_positif_nostop, iso_verif_egalite_nostop, &
37!          errmax,errmaxrel,deltalim
38       USE isotopes_mod, ONLY: ridicule,iso_eau,iso_hdo
39#endif
40#endif
41USE lmdz_cvflag
42  USE lmdz_cvthermo
43  USE lmdz_cv3param
44
45  IMPLICIT NONE
46
47  include "YOMCST2.h"
48
49!inputs:
50  INTEGER, INTENT (IN)                               :: ncum, nd, na
51  INTEGER, INTENT (IN)                               :: ntra, nloc
52  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
53  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
54  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
55  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
56  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
57  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
58  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
59  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
60  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
61  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
62  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
63  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
64  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
65  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
66  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
67#ifdef ISO
68  REAL, DIMENSION (ntraciso,nloc, na), INTENT (IN)            :: xtclw
69  REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN)            :: xt
70  REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN)            :: xtta
71#endif
72
73!outputs:
74  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
75  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
76  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
77  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax           ! Highest mixing fraction of mixed
78                                                                         ! updraughts with the sign of (h-hp)
79  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
80  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
81  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
82  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
83#ifdef ISO
84  REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT)       :: xtent
85  REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT)       :: xtelij
86#endif
87
88
89!local variables:
90  INTEGER i, j, k, il, im, jm
91  INTEGER num1, num2
92  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
93  REAL                               :: alt, delp, delm
94  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
95  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
96  REAL, DIMENSION (nloc)             :: signhpmh
97  REAL, DIMENSION (nloc)             :: Sx
98  REAL                               :: Scrit2
99  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
100  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
101  REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
102  REAL, DIMENSION (nloc, nd, nd)     :: Sij
103  REAL, DIMENSION (nloc, nd)         :: csum
104  REAL                               :: awat
105  REAL                               :: cpm        !Mixed draught heat capacity
106  REAL                               :: Tm         !Mixed draught temperature
107  LOGICAL, DIMENSION (nloc)          :: lwork
108
109  REAL amxupcrit, df, ff
110  INTEGER nstep
111
112  INTEGER,SAVE                                       :: igout=1
113!$OMP THREADPRIVATE(igout)
114
115! --   Mixing probability distribution functions
116
117  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
118
119  Qcoef1(F) = tanh(F/gammas)
120  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
121  QFF(F) = max(min(F,1.), 0.)
122  QFFf(F) = min(QFF(F), scut)
123  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
124  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
125  Qmix2(F) = -log(1.-QFFf(F))/scut
126  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
127  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
128  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
129
130  INTEGER, SAVE :: ifrst
131  DATA ifrst/0/
132!$OMP THREADPRIVATE(ifrst)
133
134#ifdef ISO
135      INTEGER ixt
136      REAL xtrti(ntraciso)
137      REAL xtres(ntraciso)
138      REAL zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc)
139      REAL qent_diag(nloc,nd,nd)
140      REAL xtent_diag(ntraciso,nloc,nd,nd)
141#endif   
142
143
144! =====================================================================
145! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
146! =====================================================================
147
148! -- Initialize mixing PDF coefficients
149  IF (ifrst==0) THEN
150    ifrst = 1
151    Qcoef1max = Qcoef1(Fmax)
152    Qcoef2max = Qcoef2(Fmax)
153!<jyg
154   PRINT*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
155            fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
156!>jyg
157
158  END IF
159
160
161! ori        do 360 i=1,ncum*nlp
162  DO j = 1, nl
163    DO i = 1, ncum
164      nent(i, j) = 0
165! in convect3, m is computed in cv3_closure
166! ori          m(i,1)=0.0
167    END DO
168  END DO
169
170#ifdef ISOVERIF
171        Qent(:,:,:)=0.0
172        xtent(:,:,:,:)=0.0
173#endif
174
175! ori      do 400 k=1,nlp
176! ori       do 390 j=1,nlp
177  DO j = 1, nl
178    DO k = 1, nl
179      DO i = 1, ncum
180        Qent(i, k, j) = rr(i, j)
181        uent(i, k, j) = u(i, j)
182        vent(i, k, j) = v(i, j)
183        elij(i, k, j) = 0.0
184        hent(i, k, j) = 0.0
185!AC!            Ment(i,k,j)=0.0
186!AC!            Sij(i,k,j)=0.0
187#ifdef ISO
188        xtelij(:,i, k, j) = 0.0
189        xtent(:,i, k, j) = xt(:,i, j)
190#endif
191      END DO
192    END DO
193  END DO
194#ifdef ISO
195        ! on initialise mieux par sécurité:
196        elij(:, :, :)=0.0
197        xtelij(:,:, :, :)=0.0
198#endif
199
200!AC!
201  Ment(1:ncum, 1:nd, 1:nd) = 0.0
202  Sij(1:ncum, 1:nd, 1:nd) = 0.0
203!AC!
204!ym
205  Sigij(1:ncum, 1:nd, 1:nd) = 0.0
206!ym
207
208!jyg!  DO k = 1, ntra
209!jyg!    DO j = 1, nd ! instead nlp
210!jyg!      DO i = 1, nd ! instead nlp
211!jyg!        DO il = 1, ncum
212!jyg!          traent(il, i, j, k) = tra(il, j, k)
213!jyg!        END DO
214!jyg!      END DO
215!jyg!    END DO
216!jyg!  END DO
217
218! =====================================================================
219! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
220! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
221! --- FRACTION (Sij)
222! =====================================================================
223
224  DO i = minorig + 1, nl
225
226    IF (ok_entrain) THEN
227      DO j = minorig, nl
228        DO il = 1, ncum
229          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
230                           .AND. (j<=inb(il))) THEN
231
232!!            rti = qnk(il) - ep(il, i)*clw(il, i)
233            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
234            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
235!jyg(from aj)<
236            IF (cvflag_ice) THEN
237! PRINT*,cvflag_ice,'cvflag_ice dans do 700'
238              IF (t(il,j)<=263.15) THEN
239                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
240                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
241              END IF
242            END IF ! IF (cvflag_ice) THEN
243!>jyg
244            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
245            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
246            dei = denom
247            IF (abs(dei)<0.01) dei = 0.01
248            Sij(il, i, j) = anum/dei
249            Sij(il, i, i) = 1.0
250            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
251            altem = altem/bf2
252            cwat = clw(il, j)*(1.-ep(il,j))
253            stemp = Sij(il, i, j)
254            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
255!jyg(from aj)<
256              IF (cvflag_ice) THEN
257                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
258                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
259              ELSE
260                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
261                denom = denom + lv(il, j)*(rr(il,i)-rti)
262              END IF
263!>jyg
264              IF (abs(denom)<0.01) denom = 0.01
265              Sij(il, i, j) = anum/denom
266              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
267              altem = altem - (bf2-1.)*cwat
268            END IF
269            IF (Sij(il,i,j)>0.0) THEN
270!!!                 Ment(il,i,j)=m(il,i)
271              Ment(il, i, j) = 1.
272              elij(il, i, j) = altem
273              elij(il, i, j) = amax1(0.0, elij(il,i,j))
274              nent(il, i) = nent(il, i) + 1
275            END IF
276
277            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
278            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
279          ELSE IF (j > i) THEN
280            IF (prt_level >= 10) THEN
281              print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
282            ENDIF
283          END IF ! new ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
284#ifdef ISO
285         IF(sij(il,i,j).gt.0.0)THEN
286               qent_diag(il,i,j)=sij(il,i,j)*rr(il,i) &
287                     +(1.-sij(il,i,j))*rti
288           ! question: pourquoi qent n'est-il pas recalculé ici?
289         endif !IF(sij(il,i,j).gt.0.0)THEN
290#endif
291        END DO ! il
292
293#ifdef ISO
294       do il=1,ncum
295         zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
296         zfice(il) = MIN(MAX(zfice(il),0.0),1.0)
297
298         IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. &
299           (j.ge.(icb(il)-1)).AND.(j.le.inb(il))) THEN
300          IF (sij(il,i,j).gt.0.0) THEN
301           do ixt=1,ntraciso
302            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
303            xtent_diag(ixt,il,i,j)=sij(il,i,j)*xt(ixt,il,i) &
304                     +(1.-sij(il,i,j))*xtrti(ixt)   
305           enddo ! do ixt=1,ntraciso
306#ifdef ISOVERIF
307          IF (iso_eau.gt.0) THEN
308          CALL iso_verif_egalite(xtta(iso_eau,il,i-1), &
309                qta(il,i-1),'cv3p_mixing 256a')
310          CALL iso_verif_egalite(xt(iso_eau,il,i), &
311                rr(il,i),'cv3p_mixing 256b')
312          CALL iso_verif_egalite(xtclw(iso_eau,il,i), &
313                clw(il,i),'cv3p_mixing 256c')
314          CALL iso_verif_egalite(xtent_diag(iso_eau,il,i,j), &
315                qent_diag(il,i,j),'cv3p_mixing 256d')   
316          endif ! if (iso_eau.gt.0) THEN
317#ifdef ISOTRAC
318          CALL iso_verif_traceur(xtrti,'cv3p_mixing 1968')
319#endif           
320#endif           
321           IF (qent_diag(il,i,j).lt.elij(il,i,j)) THEN
322             do ixt=1,ntraciso
323               xtent_diag(ixt,il,i,j)=  xtent_diag(ixt,il,i,j) &
324                / qent_diag(il,i,j)*elij(il,i,j)   
325             enddo ! do ixt=1,ntraciso
326             qent_diag(il,i,j)=elij(il,i,j)
327           endif !if (qent_diag(il,i,j).lt.elij(il,i,j)) THEN
328#ifdef ISOVERIFIDRIS
329          CALL iso_verif_noNaN(qent_diag(il,i,j),'cv3p_routines 261')
330#endif
331#ifdef ISOVERIF
332          IF (iso_eau.gt.0) THEN
333           IF (iso_verif_egalite_nostop(xtent_diag(iso_eau,il,i,j), &
334                qent_diag(il,i,j),'cv3p_mixing 258').EQ.1) THEN
335             WRITE(*,*) 'il,i,j=',il,i,j
336             WRITE(*,*) 'qent_diag(il,i,j),elij(il,i,j)=', &
337                qent_diag(il,i,j),elij(il,i,j)
338             WRITE(*,*) 'sij(il,i,j)=',sij(il,i,j)
339             WRITE(*,*) 'xt=',xt(iso_eau,il,i),rr(il,i)
340             WRITE(*,*) 'xtclw=',xtclw(iso_eau,il,i),clw(il,i)
341             WRITE(*,*) 'ep(il,i)=',ep(il,i)
342             stop
343           endif
344          endif
345          CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
346                'cv3p_routines 262')
347#endif
348         else !if (sij(il,i,j).gt.0.0) THEN
349            ! ajout le 4 mai 2012: initialiser au cas où
350            elij(il,i,j)=0.0
351            qent_diag(il,i,j)=rr(il,i)
352            do ixt=1,ntraciso
353              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
354            enddo !do ixt=1,ntraciso
355         endif  !IF(sij(il,i,j).gt.0.0)THEN
356        else  !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
357            ! ajout le 4 mai 2012: initialiser au cas où
358            elij(il,i,j)=0.0
359            qent_diag(il,i,j)=rr(il,i)
360            do ixt=1,ntraciso
361              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
362            enddo !do ixt=1,ntraciso
363        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
364#ifdef ISOVERIF
365#ifdef VERIFNEGATIF
366        CALL iso_verif_positif(qent_diag(il,i,j), &
367                'cv3p_mixing 303: qt<0')
368#endif     
369        CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
370                'cv3p_mixing 305: qt<elij')
371#endif       
372       enddo  !do il=1,ncum
373
374       !WRITE(*,*) 'cv3p_mixing 265: CALL condiso'
375       CALL condiso_liq_ice_vectall(xtent_diag(1,1,i,j), &
376                qent_diag(1,i,j),elij(1,i,j), &
377                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
378#ifdef ISOTRAC
379        CALL condiso_liq_ice_vectall_trac(xtent_diag(1,1,i,j), &
380                qent_diag(1,i,j),elij(1,i,j), &
381                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum) 
382#ifdef ISOVERIF
383        do il=1,ncum
384          CALL iso_verif_traceur(xt(1,il,i),'cv3p_mixing 1967')
385          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1969')
386        enddo !do il=1,ncum
387#endif           
388#endif   
389        do il=1,ncum
390         do ixt = 1, ntraciso
391          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
392         enddo !do ixt = 1, ntraciso
393        enddo !do il=1,ncum
394
395#ifdef ISOTRAC   
396!        WRITE(*,*) 'cv3p_mixing tmp 1987,option_traceurs=',
397!     :           option_traceurs
398        IF (option_tmin.ge.1) THEN
399        do il=1,ncum   
400        ! colorier la vapeur résiduelle selon température de
401        ! condensation, et le condensat en un tag spécifique
402          IF ((elij(il,i,j).gt.0.0).AND. &
403                (qent_diag(il,i,j).gt.0.0)) THEN
404            IF (option_traceurs.EQ.17) THEN
405             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
406                elij(il,i,j), &
407                xtent_diag(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
408                0.0,xtres, &
409                seuil_tag_tmin)
410            else !if (option_traceurs.EQ.17) THEN
411!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
412             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
413                elij(il,i,j), &
414                xtent_diag(1,il,i,j),xtelij(1,il,i,j),rs(il,j), &
415                0.0,xtres,seuil_tag_tmin)
416            endif !if (option_traceurs.EQ.17) THEN
417            do ixt=1+niso,ntraciso
418               xtent_diag(ixt,il,i,j)=xtres(ixt)
419            enddo     
420          endif !if (cond.gt.0.0) THEN
421        enddo !do il=1,ncum
422#ifdef ISOVERIF
423        do il=1,ncum
424          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1996')
425          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3p_mixing 1997')
426          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
427                'cv3p_mixing 2042')
428        enddo !do il=1,ncum
429#endif       
430        endif !if (option_tmin.ge.1) THEN
431#endif
432
433! fractionation:
434#ifdef ISOVERIF                 
435        do il=1,ncum
436        IF ((i.ge.icb(il)).AND.(i.le.inb(il)).AND. &
437           (j.ge.(icb(il)-1)).AND.(j.le.inb(il))) THEN
438        IF (sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.scut) THEN
439        IF (iso_eau.gt.0) THEN
440          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,j), &
441             qent(il,i,j),'cv3p_mixing 1889',errmax,errmaxrel)   
442          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,j), &
443             elij(il,i,j),'cv3p_mixing 1890',errmax,errmaxrel)
444        endif
445        IF (iso_hdo.gt.0) THEN
446          CALL iso_verif_aberrant_choix(xt(iso_HDO,il,i),rr(il,i), &
447                 ridicule,deltalim,'cv3p_mixing 1997')         
448          CALL iso_verif_aberrant_choix( &
449                 xtent(iso_HDO,il,i,j),qent(il,i,j), &
450                 ridicule,deltalim,'cv3p_mixing 1931')
451          CALL iso_verif_aberrant_choix( &
452                 xtelij(iso_HDO,il,i,j),elij(il,i,j), &
453                 ridicule,deltalim_snow,'cv3p_mixing 1993')
454        endif !if (iso_hdo.gt.0) THEN
455#ifdef ISOTRAC 
456!        WRITE(*,*) 'cv3p_mixing tmp 2039 il=',il
457           CALL iso_verif_traceur(xtent(1,il,i,j), &
458                        'cv3p_mixing 2031')
459           CALL iso_verif_traceur(xtelij(1,il,i,j), &
460                        'cv3p_mixing 2033')
461#endif       
462
463        endif !IF(sij(il,i,j).gt.0.0)THEN
464        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
465        enddo !do il=1,ncum
466#endif
467!        WRITE(*,*) 'cv3p_routine tmp 1984: cond=',elij(il,i,j)
468#endif
469
470      END DO ! j
471    ELSE  ! (ok_entrain)
472      DO il = 1,ncum
473        nent(il,i) = 0
474      ENDDO
475    ENDIF ! (ok_entrain)
476
477!jygdebug<
478    IF (prt_level >= 10) THEN
479      print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
480      IF (nent(igout,i) > 0) THEN
481        print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
482      ENDIF
483    ENDIF
484!>jygdebug
485
486! ***   if no air can entrain at level i assume that updraft detrains  ***
487! ***   at that level and calculate detrained air flux and properties  ***
488
489
490! @      do 170 i=icb(il),inb(il)
491
492    DO il = 1, ncum
493      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
494! @      IF(nent(il,i).EQ.0)THEN
495!!!       Ment(il,i,i)=m(il,i)
496        Ment(il, i, i) = 1.
497!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
498        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
499        uent(il, i, i) = unk(il)
500        vent(il, i, i) = vnk(il)
501        IF (fl_cor_ebil >= 2) THEN
502          hent(il, i, i) = hp(il,i)
503        ENDIF
504        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
505        Sij(il, i, i) = 0.0
506
507
508#ifdef ISO
509      do ixt = 1, ntraciso
510       xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
511       xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite
512      enddo ! do ixt = 1, ntraciso
513#ifdef ISOVERIF
514        IF (iso_eau.gt.0) THEN
515          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
516             qent(il,i,i),'cv3p_mixing 414',errmax,errmaxrel)
517        endif
518#endif
519#ifdef ISOTRAC         
520        IF (option_tmin.ge.1) THEN
521        ! colorier la vapeur résiduelle selon température de
522        ! condensation, et le condensat en un tag spécifique
523!        WRITE(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=',
524!     :            il,i,j,xtent(:,il,i,j)
525          IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN
526            IF (option_traceurs.EQ.17) THEN
527             CALL iso_recolorise_condensation(qent(il,i,i), &
528                elij(il,i,i), &
529                xt(1,il,nk(il)),xtclw(1,il,i),t(il,i),ep(il,i), &
530                xtres, &
531                seuil_tag_tmin)
532            else !if (option_traceurs.EQ.17) THEN
533             CALL iso_recolorise_condensation(qent(il,i,i), &
534                elij(il,i,i), &
535                xt(1,il,nk(il)),xtclw(1,il,i),rs(il,i),ep(il,i), &
536                xtres, &
537                seuil_tag_tmin)
538            endif !if (option_traceurs.EQ.17) THEN
539            do ixt=1+niso,ntraciso
540              xtent(ixt,il,i,i)=xtres(ixt)
541            enddo
542#ifdef ISOVERIF           
543            do ixt=1,niso
544            CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
545                'cv3p_mixing 2102',errmax,errmaxrel)
546            CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
547                'cv3p_mixing 2154')
548            enddo
549#endif           
550          endif !if (cond.gt.0.0) THEN
551#ifdef ISOVERIF         
552          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
553                qent(il,i,i),'cv3p_mixing 2103',errmax,errmaxrel)
554          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2095')
555          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2096')
556#endif       
557        endif !if (option_tmin.ge.1) THEN
558#endif
559! END IF ISOTRAC
560#endif
561! END IF ISO
562
563      END IF ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
564    END DO ! il
565  END DO ! i = minorig + 1, nl
566
567!jyg!  DO j = 1, ntra
568!jyg!    DO i = minorig + 1, nl
569!jyg!      DO il = 1, ncum
570!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
571!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
572!jyg!        END IF
573!jyg!      END DO
574!jyg!    END DO
575!jyg!  END DO
576
577  DO j = minorig, nl
578    DO i = minorig, nl
579      DO il = 1, ncum
580        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
581            (i>=icb(il)) .AND. (i<=inb(il))) THEN
582          Sigij(il, i, j) = Sij(il, i, j)
583        END IF
584      END DO
585    END DO
586  END DO
587! @      enddo
588
589! @170   continue
590
591! =====================================================================
592! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
593! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
594! =====================================================================
595
596  CALL zilch(csum, nloc*nd)
597
598  DO il = 1, ncum
599    lwork(il) = .FALSE.
600  END DO
601
602! ---------------------------------------------------------------
603  DO i = minorig + 1, nl      !Loop on origin level "i"
604! ---------------------------------------------------------------
605
606    num1 = 0
607    DO il = 1, ncum
608      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
609    END DO
610    IF (num1<=0) GO TO 789
611
612
613!JYG1    Find maximum of SIJ for J>I, if any.
614
615    Sx(:) = 0.
616
617    DO il = 1, ncum
618      IF (i>=icb(il) .AND. i<=inb(il)) THEN
619        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
620        Sbef(il) = max(0., signhpmh(il))
621      END IF
622    END DO
623
624    DO j = i + 1, nl
625      DO il = 1, ncum
626        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
627          IF (Sbef(il)<Sij(il,i,j)) THEN
628            Sx(il) = max(Sij(il,i,j), Sx(il))
629          END IF
630          Sbef(il) = Sij(il, i, j)
631        END IF
632      END DO
633    END DO
634
635
636    DO il = 1, ncum
637      IF (i>=icb(il) .AND. i<=inb(il)) THEN
638        lwork(il) = (nent(il,i)/=0)
639!!        rti = qnk(il) - ep(il, i)*clw(il, i)
640        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
641!jyg<
642        IF (cvflag_ice) THEN
643
644          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
645                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
646          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
647                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
648        ELSE
649
650          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
651                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
652          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
653                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
654        END IF
655!>jyg
656        IF (abs(denom)<0.01) denom = 0.01
657        Scrit(il) = min(anum/denom, 1.)
658        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
659
660!JYG1    Find new critical value Scrit2
661!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
662!                     Sij < Scrit2  => mixed draught will detrain at J>I
663
664        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
665                 Scrit(il)*max(0., signhpmh(il))
666
667        Scrit(il) = Scrit2
668
669!JYG    Correction pour la nouvelle logique; la correction pour ALT
670! est un peu au hazard
671        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
672        IF (alt<=0.0) Scrit(il) = 1.0
673
674        smax(il) = 0.0
675        ASij(il) = 0.0
676        sup(il) = 0.      ! upper S-value reached by descending draughts
677      END IF
678    END DO
679
680! ---------------------------------------------------------------
681    DO j = minorig, nl         !Loop on destination level "j"
682! ---------------------------------------------------------------
683
684      num2 = 0
685      DO il = 1, ncum
686        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
687            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
688            lwork(il)) num2 = num2 + 1
689      END DO
690      IF (num2<=0) GO TO 175
691
692! -----------------------------------------------
693      IF (j>i) THEN
694! -----------------------------------------------
695        DO il = 1, ncum
696          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
697              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
698              lwork(il)) THEN
699            IF (Sij(il,i,j)>0.0) THEN
700              Smid(il) = min(Sij(il,i,j), Scrit(il))
701              Sjmax(il) = Smid(il)
702              Sjmin(il) = Smid(il)
703              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
704                smin(il) = Smid(il)
705                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
706                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
707                Sjmin(il) = min(Sjmin(il), Scrit(il))
708                Sbef(il) = Sij(il, i, j)
709              END IF
710            END IF
711          END IF
712        END DO
713! -----------------------------------------------
714      ELSE IF (j==i) THEN
715! -----------------------------------------------
716        DO il = 1, ncum
717          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
718              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
719              lwork(il)) THEN
720            IF (Sij(il,i,j)>0.0) THEN
721              Smid(il) = 1.
722              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
723                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
724              Sjmin(il) = max(Sjmin(il), sup(il))
725              Sjmax(il) = 1.
726
727! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
728              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
729
730              smin(il) = 1.
731              Sbef(il) = max(0., signhpmh(il))
732              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
733            END IF
734          END IF
735        END DO
736! -----------------------------------------------
737      ELSE IF (j<i) THEN
738! -----------------------------------------------
739        DO il = 1, ncum
740          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
741              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
742              lwork(il)) THEN
743            IF (Sij(il,i,j)>0.0) THEN
744              Smid(il) = max(Sij(il,i,j), Scrit(il))
745              Sjmax(il) = Smid(il)
746              Sjmin(il) = Smid(il)
747              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
748                smax(il) = Smid(il)
749                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
750                Sjmax(il) = max(Sjmax(il), Scrit(il))
751                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
752                Sjmin(il) = max(Sjmin(il), Scrit(il))
753                Sbef(il) = Sij(il, i, j)
754              END IF
755              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
756                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
757            END IF
758          END IF
759        END DO
760! -----------------------------------------------
761      END IF
762! -----------------------------------------------
763
764
765      DO il = 1, ncum
766        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
767            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
768            lwork(il)) THEN
769          IF (Sij(il,i,j)>0.0) THEN
770!!            rti = qnk(il) - ep(il, i)*clw(il, i)
771            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
772            Qmixmax(il) = Qmix(Sjmax(il))
773            Qmixmin(il) = Qmix(Sjmin(il))
774            Rmixmax(il) = Rmix(Sjmax(il))
775            Rmixmin(il) = Rmix(Sjmin(il))
776            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
777            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
778
779            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
780
781! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
782            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
783              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
784            ELSE
785              Sigij(il, i, j) = 0.
786            END IF
787
788! --    Compute Qent, uent, vent according to the true mixing fraction
789            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
790            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
791            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
792
793! --     Compute liquid water static energy of mixed draughts
794!    IF (j .GT. i) THEN
795!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
796!      awat=amax1(awat,0.0)
797!    ELSE
798!      awat = 0.
799!    ENDIF
800!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
801!    :         + Sigij(il,i,j)*H(il,i)
802!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
803!IM 301008 beg
804            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
805
806!jyg<
807!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
808!            elij(il, i, j) = elij(il, i, j) + &
809!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
810!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
811!            elij(il, i, j) = elij(il, i, j) / &
812!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
813!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
814
815!       Computation of condensate amount Elij, taking into account the ice fraction frac
816!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
817!                 should be corrected.
818
819!  Heat capacity of mixed draught
820    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
821
822    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
823            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
824            elij(il, i, j) = elij(il, i, j) + &
825                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
826                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
827            elij(il, i, j) = elij(il, i, j) / &
828                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
829                              (cpm*rrv*t(il,j)*t(il,j)))
830    ELSE
831            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
832            elij(il, i, j) = elij(il, i, j) + &
833                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
834                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
835            elij(il, i, j) = elij(il, i, j) / &
836                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
837                              (cpm*rrv*t(il,j)*t(il,j)))
838    ENDIF
839!>jyg
840            elij(il, i, j) = max(elij(il,i,j), 0.)
841
842            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
843
844            IF (j>i) THEN
845              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
846              awat = amax1(awat, 0.0)
847            ELSE
848              awat = 0.
849            END IF
850
851! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
852! :         t(il,j))
853
854!jyg<
855!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
856! Mixed draught temperature at level j
857    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
858          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))
859          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
860    ELSE
861          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))
862          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
863    ENDIF
864!>jyg
865
866!IM 301008 end
867
868! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
869! :               i,j,hent(il,i,j),Sigij(il,i,j)
870
871! --      ASij is the integral of P(F) over the relevant F interval
872            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
873                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
874
875          END IF
876        END IF
877      END DO ! il
878
879
880             
881#ifdef ISO
882      do il=1,ncum
883          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
884          zfice(il) = MIN(MAX(zfice(il),0.0),1.0) 
885       IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
886           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
887           .AND. lwork(il) ) THEN
888        IF(sij(il,i,j).gt.0.0)THEN
889          do ixt=1,ntraciso           
890            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
891            xtent(ixt,il,i,j) = (1.-Sigij(il,i,j))*xtrti(ixt) &
892                    + Sigij(il,i,j)*xt(ixt,il,i)
893          enddo !do ixt = 1, ntraciso     
894#ifdef ISOVERIF
895#ifdef ISOTRAC         
896          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
897          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
898#endif         
899          IF (iso_eau.gt.0) THEN
900            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
901                qent(il,i,j),'cv3p_routines 718')   
902          endif
903#endif             
904         endif  !IF(sij(il,i,j).gt.0.0)THEN
905        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
906       enddo  !do il=1,ncum
907   
908#ifdef ISOVERIFIDRIS
909        do il=1,ncum       
910         CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 771')
911        enddo !do il=1,ncum   
912#endif     
913#ifdef ISOVERIF
914        IF (iso_eau.gt.0) THEN
915         do il=1,ncum   
916            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
917                qent(il,i,j),'cv3p_routines 744')   
918         enddo !do il=1,ncum   
919        endif   
920        do il=1,ncum       
921         CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
922                'cv3p_routines 749')
923        enddo !do il=1,ncum   
924        !WRITE(*,*) 'cv3p_mixing 788: CALL condiso_liq_ice_vectall'
925#endif     
926     
927       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
928                elij(1,i,j), &
929                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
930#ifdef ISOTRAC
931        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
932                elij(1,i,j), &
933                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
934#endif   
935        do il=1,ncum
936        IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
937           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
938           .AND. lwork(il) ) THEN
939        IF(sij(il,i,j).gt.0.0)THEN
940         do ixt = 1, ntraciso
941          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
942         enddo !do ixt = 1, ntraciso
943        endif !IF(sij(il,i,j).gt.0.0)THEN
944        endif !if ( i.ge.icb(il) .AND. i.le.inb(il) .AND.
945        enddo !do il=1,ncum
946
947#ifdef ISOVERIF 
948        IF (iso_eau.gt.0) THEN
949            do il=1,ncum
950            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
951                qent(il,i,j),'cv3p_routines 744')
952            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
953                elij(il,i,j),'cv3p_routines 746')
954            enddo !do il=1,ncum   
955        endif
956#endif 
957!#ifdef ISOVERIF 
958
959#ifdef ISOTRAC   
960!        WRITE(*,*) 'cv3_routines tmp 1987,option_traceurs=',
961!     :           option_traceurs
962        IF (option_tmin.ge.1) THEN
963        do il=1,ncum   
964!        WRITE(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
965!     :           'tcond(il),rs(il,j)=',
966!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
967        ! colorier la vapeur résiduelle selon température de
968        ! condensation, et le condensat en un tag spécifique
969          IF ((elij(il,i,j).gt.0.0).AND.(qent(il,i,j).gt.0.0)) THEN
970            IF (option_traceurs.EQ.17) THEN
971             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
972                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
973                0.0,xtres, &
974                seuil_tag_tmin)
975            else !if (option_traceurs.EQ.17) THEN
976!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
977             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
978                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
979                seuil_tag_tmin)
980            endif !if (option_traceurs.EQ.17) THEN
981            do ixt=1+niso,ntraciso
982               xtent(ixt,il,i,j)=xtres(ixt)
983            enddo     
984          endif !if (cond.gt.0.0) THEN
985        enddo !do il=1,ncum
986#ifdef ISOVERIF
987        do il=1,ncum
988          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
989          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
990          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
991                'cv3_routines 2042')
992        enddo !do il=1,ncum
993#endif       
994        endif !if (option_tmin.ge.1) THEN
995#endif
996! END IF ISOTRAC
997
998#endif
999!!#ifdef ISO
1000
1001!jyg!      DO k = 1, ntra
1002!jyg!        DO il = 1, ncum
1003!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1004!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1005!jyg!              lwork(il)) THEN
1006!jyg!            IF (Sij(il,i,j)>0.0) THEN
1007!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
1008!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
1009!jyg!            END IF
1010!jyg!          END IF
1011!jyg!        END DO
1012!jyg!      END DO
1013
1014! --    If I=J (detrainement and entrainement at the same level), then only the
1015! --    adiabatic ascent part of the mixture is considered
1016      IF (i==j) THEN
1017        DO il = 1, ncum
1018          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
1019              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
1020              lwork(il)) THEN
1021            IF (Sij(il,i,j)>0.0) THEN
1022!!              rti = qnk(il) - ep(il, i)*clw(il, i)
1023              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
1024!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
1025              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
1026                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
1027              Qent(il, i, i) = rti
1028              uent(il, i, i) = unk(il)
1029              vent(il, i, i) = vnk(il)
1030              hent(il, i, i) = hp(il, i)
1031              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1032              Sigij(il, i, i) = 0.
1033            END IF
1034          END IF
1035        END DO
1036
1037
1038#ifdef ISO
1039      IF (I .EQ. J) THEN
1040      do il=1,ncum
1041          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
1042          zfice(il) = MIN(MAX(zfice(il),0.0),1.0) 
1043       IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
1044           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
1045           .AND. lwork(il) ) THEN
1046        IF(sij(il,i,j).gt.0.0)THEN
1047          do ixt=1,ntraciso           
1048            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
1049            xtent(ixt,il,i,j) = xtrti(ixt)
1050          enddo !do ixt = 1, ntraciso
1051#ifdef ISOVERIFIDRIS
1052          CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 923')
1053#endif
1054#ifdef ISOVERIF 
1055#ifdef ISOTRAC         
1056          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
1057          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
1058#endif         
1059          IF (iso_eau.gt.0) THEN
1060            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
1061                qent(il,i,j),'cv3p_routines 718')   
1062          endif         
1063          CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
1064                'cv3p_routines 887')
1065#endif
1066         endif  !IF(sij(il,i,j).gt.0.0)THEN
1067        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
1068       enddo  !do il=1,ncum
1069     
1070       !WRITE(*,*) 'cv3p_mixing 808: CALL condiso'
1071       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
1072                elij(1,i,j), &
1073                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1074
1075       do il=1,ncum
1076         do ixt = 1, ntraciso
1077          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
1078         enddo !do ixt = 1, ntraciso
1079#ifdef ISOVERIF
1080         IF (iso_eau.gt.0) THEN
1081           CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1082                elij(il,i,j),'cv3p_routines 1074')
1083         endif ! if (iso_eau.gt.0) THEN
1084#endif
1085        enddo !do il=1,ncum
1086       
1087#ifdef ISOTRAC
1088        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
1089                elij(1,i,j), &
1090                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1091#ifdef ISOVERIF
1092        do il=1,ncum
1093          CALL iso_verif_traceur(xt(1,il,i),'cv3p_routines 1967')
1094          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 1969')
1095        enddo !do il=1,ncum
1096#endif   
1097       
1098       IF (option_tmin.ge.1) THEN
1099        do il=1,ncum   
1100!        WRITE(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
1101!     :           'tcond(il),rs(il,j)=',
1102!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
1103        ! colorier la vapeur résiduelle selon température de
1104        ! condensation, et le condensat en un tag spécifique
1105          IF ((elij(il,i,j).gt.0.0).AND.(qent(il,i,j).gt.0.0)) THEN
1106            IF (option_traceurs.EQ.17) THEN
1107             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1108                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
1109                0.0,xtres, &
1110                seuil_tag_tmin)
1111            else !if (option_traceurs.EQ.17) THEN
1112!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
1113             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1114                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
1115                seuil_tag_tmin)
1116            endif !if (option_traceurs.EQ.17) THEN
1117            do ixt=1+niso,ntraciso
1118               xtent(ixt,il,i,j)=xtres(ixt)
1119            enddo     
1120          endif !if (cond.gt.0.0) THEN
1121        enddo !do il=1,ncum
1122#ifdef ISOVERIF
1123        do il=1,ncum
1124          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
1125          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
1126          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1127                'cv3_routines 2042')
1128        enddo !do il=1,ncum
1129#endif       
1130        endif !if (option_tmin.ge.1) THEN
1131!#ifdef ISOVERIF       
1132#endif
1133!#ifdef ISOTRAC       
1134#ifdef ISOVERIF       
1135        IF (iso_eau.gt.0) THEN
1136            do il=1,ncum
1137            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
1138                qent(il,i,j),'cv3p_routines 883')
1139            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1140                elij(il,i,j),'cv3p_routines 885') 
1141            enddo !do il=1,ncum 
1142        endif !if (iso_eau.gt.0) THEN
1143#endif 
1144!#ifdef ISOVERIF           
1145        endif !IF (I .EQ. J) THEN
1146#endif 
1147!#ifdef ISO
1148
1149!jyg!        DO k = 1, ntra
1150!jyg!          DO il = 1, ncum
1151!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1152!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1153!jyg!                lwork(il)) THEN
1154!jyg!              IF (Sij(il,i,j)>0.0) THEN
1155!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
1156!jyg!              END IF
1157!jyg!            END IF
1158!jyg!          END DO
1159!jyg!        END DO
1160
1161      END IF
1162
1163! ---------------------------------------------------------------
1164175 END DO        ! End loop on destination level "j"
1165! ---------------------------------------------------------------
1166
1167    DO il = 1, ncum
1168      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1169        ASij(il) = amax1(1.0E-16, ASij(il))
1170!jyg+lluis<
1171!!        ASij(il) = 1.0/ASij(il)
1172        ASij_inv(il) = 1.0/ASij(il)
1173!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
1174        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
1175!>jyg+lluis
1176        csum(il, i) = 0.0
1177      END IF
1178    END DO
1179
1180    DO j = minorig, nl
1181      DO il = 1, ncum
1182        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1183            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1184!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
1185          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
1186        END IF
1187      END DO
1188    END DO
1189
1190    DO j = minorig, nl
1191      DO il = 1, ncum
1192        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1193            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1194          csum(il, i) = csum(il, i) + Ment(il, i, j)
1195        END IF
1196      END DO
1197    END DO
1198
1199    DO il = 1, ncum
1200      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1201! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
1202        nent(il, i) = 0
1203! cc        Ment(il,i,i)=m(il,i)
1204        Ment(il, i, i) = 1.
1205!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1206        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
1207        uent(il, i, i) = unk(il)
1208        vent(il, i, i) = vnk(il)
1209        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1210
1211#ifdef ISO
1212      do ixt = 1, ntraciso
1213        xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
1214        xtelij(ixt,il,i,i)=xtclw(ixt,il,i)*(1.-ep(il,i))
1215      enddo
1216#ifdef ISOVERIF
1217      IF (iso_eau.gt.0) THEN
1218        CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1219                elij(il,i,j),'cv3p_routines 1204')
1220      endif ! if (iso_eau.gt.0) THEN
1221#endif
1222#endif
1223
1224#ifdef ISOTRAC         
1225        IF (option_tmin.ge.1) THEN
1226        ! colorier la vapeur résiduelle selon température de
1227        ! condensation, et le condensat en un tag spécifique
1228!        WRITE(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=',
1229!     :            il,i,j,xtent(:,il,i,j)
1230          IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN
1231            IF (option_traceurs.EQ.17) THEN
1232              CALL iso_recolorise_condensation(qent(il,i,i), &
1233                elij(il,i,i), &
1234                xt(1,il,1),xtclw(1,il,i),t(il,i),ep(il,i), &
1235                xtres, &
1236                seuil_tag_tmin)
1237            else !if (option_traceurs.EQ.17) THEN
1238              CALL iso_recolorise_condensation(qent(il,i,i), &
1239                elij(il,i,i), &
1240                xt(1,il,1),xtclw(1,il,i),rs(il,i),ep(il,i), &
1241                xtres, &
1242                seuil_tag_tmin)
1243            endif ! if (option_traceurs.EQ.17) THEN
1244            do ixt=1+niso,ntraciso
1245              xtent(ixt,il,i,i)=xtres(ixt)
1246            enddo 
1247#ifdef ISOVERIF               
1248            do ixt=1,niso
1249              CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
1250                'cv3p_mixing 2318',errmax,errmaxrel)
1251              CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1252                'cv3p_mixing 2383')
1253            enddo
1254#endif               
1255          endif !if (cond.gt.0.0) THEN
1256#ifdef ISOVERIF         
1257          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
1258                qent(il,i,i),'cv3p_mixing 2321',errmax,errmaxrel)
1259          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), &
1260                elij(il,i,i),'cv3p_mixing 2321b',errmax,errmaxrel)
1261          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2322')
1262          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2323')
1263#endif       
1264        endif !if (option_tmin.ge.1) THEN
1265#endif
1266        IF (fl_cor_ebil .GE. 2) THEN
1267          hent(il, i, i) = hp(il,i)
1268          Sigij(il, i, i) = 0.0
1269        ELSE
1270          Sij(il, i, i) = 0.0
1271        ENDIF
1272      END IF
1273    END DO ! il
1274
1275!jyg!    DO j = 1, ntra
1276!jyg!      DO il = 1, ncum
1277!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1278!jyg!! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
1279!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
1280!jyg!        END IF
1281!jyg!      END DO
1282!jyg!    END DO
1283
1284! ---------------------------------------------------------------
1285789 END DO              ! End loop on origin level "i"
1286! ---------------------------------------------------------------
1287
1288#ifdef ISO
1289#ifdef ISOVERIF
1290!       WRITE(*,*) 'cv3p_mixing 2540: ', &
1291!        'verif finale en sortant de cv3p_mixing'
1292!       WRITE(*,*) 'qent,xtent(1,1,1)=',qent(1,1,1),xtent(iso_eau,1,1,1)
1293       do im = 1, nd
1294       do jm = 1, nd
1295        do il = 1, ncum
1296          IF (iso_eau.gt.0) THEN
1297            CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1298              elij(il,im,jm),'cv3p_mixing 2110',errmax,errmaxrel)
1299            CALL iso_verif_egalite_choix(xtent(iso_eau,il,im,jm), &
1300              qent(il,im,jm),'cv3p_mixing 2112',errmax,errmaxrel)
1301          endif !if (iso_eau>0) THEN
1302#ifdef ISOTRAC
1303        CALL iso_verif_traceur_justmass(xtelij(1,il,im,jm), &
1304                       'cv3p_routine 2250')
1305#endif           
1306        enddo !do il = 1, nloc
1307       enddo !do jm = 1, klev
1308       enddo !do im = 1, klev
1309#endif
1310#endif 
1311
1312#ifdef ISO
1313#ifdef ISOTRAC
1314        ! seulement à la fin on taggue le condensat
1315        IF (option_cond.ge.1) THEN
1316         do im = 1, nd
1317         do jm = 1, nd
1318         do il = 1, ncum   
1319           ! colorier le condensat en un tag spécifique
1320           do ixt=niso+1,ntraciso
1321             IF (index_zone(ixt).EQ.izone_cond) THEN
1322                xtelij(ixt,il,im,jm)=xtelij(index_iso(ixt),il,im,jm)
1323             else !if (index_zone(ixt).EQ.izone_cond) THEN
1324                xtelij(ixt,il,im,jm)=0.0
1325             endif !if (index_zone(ixt).EQ.izone_cond) THEN
1326           enddo !do ixt=1,ntraciso     
1327#ifdef ISOVERIF
1328        CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1329                elij(il,im,jm),'cv3p_mixing 2408',errmax,errmaxrel)
1330        CALL iso_verif_traceur(xtelij(1,il,im,jm), &
1331               'cv3p_mixing 358')
1332#endif     
1333         enddo !do il = 1, ncum   
1334         enddo !do jm = 1, nd
1335         enddo !do im = 1, nd
1336        ! 17 juin 2020: on supprime ci-dessous car normalement xtclw a déjà été
1337        ! recolorié ds cv3_undilute
1338!         do im = 1, nd
1339!         do il = 1, ncum   
1340!           ! colorier le condensat en un tag spécifique
1341!           do ixt=niso+1,ntraciso
1342!             if (index_zone(ixt).EQ.izone_cond) THEN
1343!                xtclw(ixt,il,im)=xtclw(index_iso(ixt),il,im)
1344!             else !if (index_zone(ixt).EQ.izone_cond) THEN
1345!                xtclw(ixt,il,im)=0.0
1346!             endif !if (index_zone(ixt).EQ.izone_cond) THEN
1347!           enddo !do ixt=1,ntraciso     
1348!#ifdef ISOVERIF
1349!        CALL iso_verif_egalite_choix(xtclw(iso_eau,il,im), &
1350!                clw(il,im),'cv3p_mixing 2427',errmax,errmaxrel)
1351!        CALL iso_verif_traceur(xtclw(1,il,im), &
1352!               'cv3p_mixing 358')
1353!        if (iso_verif_positif_nostop(xtclw(itZonIso( &
1354!                izone_cond,iso_eau),i,k)-xtclw(iso_eau,i,k) &
1355!                ,'cv3p_mixing 909').EQ.1) THEN
1356!               WRITE(*,*) 'i,k=',i,k
1357!               WRITE(*,*) 'xtclw=',xtclw(:,i,k)
1358!               WRITE(*,*) 'niso,ntraciso,index_zone,izone_cond=', &
1359!                  niso,ntraciso,index_zone,izone_cond     
1360!               stop
1361!         endif !if (iso_verif_positif_nostop(xtclw(itZonIso(
1362!#endif             
1363!         enddo !do il = 1, ncum   
1364!         enddo !do im = 1, nd
1365!         WRITE(*,*) 'xtclw(:,1,2)=',xtclw(:,1,2)
1366        endif !if (option_tmin.EQ.1) THEN
1367#endif
1368#endif         
1369
1370
1371
1372END SUBROUTINE cv3p_mixing
1373
Note: See TracBrowser for help on using the repository browser.