source: LMDZ6/trunk/libf/phylmdiso/cv3p_mixing.F90 @ 5423

Last change on this file since 5423 was 5304, checked in by abarral, 2 months ago

Turn YOMCST2.h.h into module

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