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

Last change on this file since 5367 was 5160, checked in by abarral, 4 months ago

Put .h into modules

  • Property svn:keywords set to Id
File size: 50.6 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  USE lmdz_yomcst2
45
46  IMPLICIT NONE
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#ifdef ISOVERIFIDRIS
328          CALL iso_verif_noNaN(qent_diag(il,i,j),'cv3p_routines 261')
329#endif
330#ifdef ISOVERIF
331          IF (iso_eau.gt.0) THEN
332           IF (iso_verif_egalite_nostop(xtent_diag(iso_eau,il,i,j), &
333                qent_diag(il,i,j),'cv3p_mixing 258').EQ.1) THEN
334             WRITE(*,*) 'il,i,j=',il,i,j
335             WRITE(*,*) 'qent_diag(il,i,j),elij(il,i,j)=', &
336                qent_diag(il,i,j),elij(il,i,j)
337             WRITE(*,*) 'sij(il,i,j)=',sij(il,i,j)
338             WRITE(*,*) 'xt=',xt(iso_eau,il,i),rr(il,i)
339             WRITE(*,*) 'xtclw=',xtclw(iso_eau,il,i),clw(il,i)
340             WRITE(*,*) 'ep(il,i)=',ep(il,i)
341             stop
342           endif
343          endif
344          CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
345                'cv3p_routines 262')
346#endif
347         else !if (sij(il,i,j).gt.0.0) THEN
348            ! ajout le 4 mai 2012: initialiser au cas où
349            elij(il,i,j)=0.0
350            qent_diag(il,i,j)=rr(il,i)
351            DO ixt=1,ntraciso
352              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
353            enddo !do ixt=1,ntraciso
354         endif  !IF(sij(il,i,j).gt.0.0)THEN
355        else  !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
356            ! ajout le 4 mai 2012: initialiser au cas où
357            elij(il,i,j)=0.0
358            qent_diag(il,i,j)=rr(il,i)
359            DO ixt=1,ntraciso
360              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
361            enddo !do ixt=1,ntraciso
362        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
363#ifdef ISOVERIF
364#ifdef VERIFNEGATIF
365        CALL iso_verif_positif(qent_diag(il,i,j), &
366                'cv3p_mixing 303: qt<0')
367#endif     
368        CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
369                'cv3p_mixing 305: qt<elij')
370#endif       
371       enddo  !do il=1,ncum
372
373       !WRITE(*,*) 'cv3p_mixing 265: CALL condiso'
374       CALL condiso_liq_ice_vectall(xtent_diag(1,1,i,j), &
375                qent_diag(1,i,j),elij(1,i,j), &
376                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
377#ifdef ISOTRAC
378        CALL condiso_liq_ice_vectall_trac(xtent_diag(1,1,i,j), &
379                qent_diag(1,i,j),elij(1,i,j), &
380                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum) 
381#ifdef ISOVERIF
382        DO il=1,ncum
383          CALL iso_verif_traceur(xt(1,il,i),'cv3p_mixing 1967')
384          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1969')
385        enddo !do il=1,ncum
386#endif           
387#endif   
388        DO il=1,ncum
389         DO ixt = 1, ntraciso
390          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
391         enddo !do ixt = 1, ntraciso
392        enddo !do il=1,ncum
393
394#ifdef ISOTRAC   
395!        WRITE(*,*) 'cv3p_mixing tmp 1987,option_traceurs=',
396!     :           option_traceurs
397        IF (option_tmin.ge.1) THEN
398        DO il=1,ncum
399        ! colorier la vapeur résiduelle selon température de
400        ! condensation, et le condensat en un tag spécifique
401          IF ((elij(il,i,j).gt.0.0).AND. &
402                (qent_diag(il,i,j).gt.0.0)) THEN
403            IF (option_traceurs.EQ.17) THEN
404             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
405                elij(il,i,j), &
406                xtent_diag(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
407                0.0,xtres, &
408                seuil_tag_tmin)
409            else !if (option_traceurs.EQ.17) THEN
410!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
411             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
412                elij(il,i,j), &
413                xtent_diag(1,il,i,j),xtelij(1,il,i,j),rs(il,j), &
414                0.0,xtres,seuil_tag_tmin)
415            endif !if (option_traceurs.EQ.17) THEN
416            DO ixt=1+niso,ntraciso
417               xtent_diag(ixt,il,i,j)=xtres(ixt)
418            enddo     
419          endif !if (cond.gt.0.0) THEN
420        enddo !do il=1,ncum
421#ifdef ISOVERIF
422        DO il=1,ncum
423          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1996')
424          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3p_mixing 1997')
425          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
426                'cv3p_mixing 2042')
427        enddo !do il=1,ncum
428#endif       
429        endif !if (option_tmin.ge.1) THEN
430#endif
431
432! fractionation:
433#ifdef ISOVERIF                 
434        DO il=1,ncum
435        IF ((i.ge.icb(il)).AND.(i.le.inb(il)).AND. &
436           (j.ge.(icb(il)-1)).AND.(j.le.inb(il))) THEN
437        IF (sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.scut) THEN
438        IF (iso_eau.gt.0) THEN
439          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,j), &
440             qent(il,i,j),'cv3p_mixing 1889',errmax,errmaxrel)   
441          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,j), &
442             elij(il,i,j),'cv3p_mixing 1890',errmax,errmaxrel)
443        endif
444        IF (iso_hdo.gt.0) THEN
445          CALL iso_verif_aberrant_choix(xt(iso_HDO,il,i),rr(il,i), &
446                 ridicule,deltalim,'cv3p_mixing 1997')         
447          CALL iso_verif_aberrant_choix( &
448                 xtent(iso_HDO,il,i,j),qent(il,i,j), &
449                 ridicule,deltalim,'cv3p_mixing 1931')
450          CALL iso_verif_aberrant_choix( &
451                 xtelij(iso_HDO,il,i,j),elij(il,i,j), &
452                 ridicule,deltalim_snow,'cv3p_mixing 1993')
453        endif !if (iso_hdo.gt.0) THEN
454#ifdef ISOTRAC 
455!        WRITE(*,*) 'cv3p_mixing tmp 2039 il=',il
456           CALL iso_verif_traceur(xtent(1,il,i,j), &
457                        'cv3p_mixing 2031')
458           CALL iso_verif_traceur(xtelij(1,il,i,j), &
459                        'cv3p_mixing 2033')
460#endif       
461
462        endif !IF(sij(il,i,j).gt.0.0)THEN
463        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
464        enddo !do il=1,ncum
465#endif
466!        WRITE(*,*) 'cv3p_routine tmp 1984: cond=',elij(il,i,j)
467#endif
468
469      END DO ! j
470    ELSE  ! (ok_entrain)
471      DO il = 1,ncum
472        nent(il,i) = 0
473      ENDDO
474    ENDIF ! (ok_entrain)
475
476!jygdebug<
477    IF (prt_level >= 10) THEN
478      PRINT *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
479      IF (nent(igout,i) > 0) THEN
480        PRINT *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
481      ENDIF
482    ENDIF
483!>jygdebug
484
485! ***   if no air can entrain at level i assume that updraft detrains  ***
486! ***   at that level and calculate detrained air flux and properties  ***
487
488
489! @      do 170 i=icb(il),inb(il)
490
491    DO il = 1, ncum
492      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
493! @      IF(nent(il,i).EQ.0)THEN
494!!!       Ment(il,i,i)=m(il,i)
495        Ment(il, i, i) = 1.
496!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
497        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
498        uent(il, i, i) = unk(il)
499        vent(il, i, i) = vnk(il)
500        IF (fl_cor_ebil >= 2) THEN
501          hent(il, i, i) = hp(il,i)
502        ENDIF
503        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
504        Sij(il, i, i) = 0.0
505
506
507#ifdef ISO
508      DO ixt = 1, ntraciso
509       xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
510       xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite
511      enddo ! do ixt = 1, ntraciso
512#ifdef ISOVERIF
513        IF (iso_eau.gt.0) THEN
514          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
515             qent(il,i,i),'cv3p_mixing 414',errmax,errmaxrel)
516        endif
517#endif
518#ifdef ISOTRAC         
519        IF (option_tmin.ge.1) THEN
520        ! colorier la vapeur résiduelle selon température de
521        ! condensation, et le condensat en un tag spécifique
522!        WRITE(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=',
523!     :            il,i,j,xtent(:,il,i,j)
524          IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN
525            IF (option_traceurs.EQ.17) THEN
526             CALL iso_recolorise_condensation(qent(il,i,i), &
527                elij(il,i,i), &
528                xt(1,il,nk(il)),xtclw(1,il,i),t(il,i),ep(il,i), &
529                xtres, &
530                seuil_tag_tmin)
531            else !if (option_traceurs.EQ.17) THEN
532             CALL iso_recolorise_condensation(qent(il,i,i), &
533                elij(il,i,i), &
534                xt(1,il,nk(il)),xtclw(1,il,i),rs(il,i),ep(il,i), &
535                xtres, &
536                seuil_tag_tmin)
537            endif !if (option_traceurs.EQ.17) THEN
538            DO ixt=1+niso,ntraciso
539              xtent(ixt,il,i,i)=xtres(ixt)
540            enddo
541#ifdef ISOVERIF           
542            DO ixt=1,niso
543            CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
544                'cv3p_mixing 2102',errmax,errmaxrel)
545            CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
546                'cv3p_mixing 2154')
547            enddo
548#endif           
549          endif !if (cond.gt.0.0) THEN
550#ifdef ISOVERIF         
551          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
552                qent(il,i,i),'cv3p_mixing 2103',errmax,errmaxrel)
553          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2095')
554          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2096')
555#endif       
556        endif !if (option_tmin.ge.1) THEN
557#endif
558! END IF ISOTRAC
559#endif
560! END IF ISO
561
562      END IF ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
563    END DO ! il
564  END DO ! i = minorig + 1, nl
565
566!jyg!  DO j = 1, ntra
567!jyg!    DO i = minorig + 1, nl
568!jyg!      DO il = 1, ncum
569!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
570!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
571!jyg!        END IF
572!jyg!      END DO
573!jyg!    END DO
574!jyg!  END DO
575
576  DO j = minorig, nl
577    DO i = minorig, nl
578      DO il = 1, ncum
579        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
580            (i>=icb(il)) .AND. (i<=inb(il))) THEN
581          Sigij(il, i, j) = Sij(il, i, j)
582        END IF
583      END DO
584    END DO
585  END DO
586! @      enddo
587
588! @170   continue
589
590! =====================================================================
591! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
592! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
593! =====================================================================
594
595  CALL zilch(csum, nloc*nd)
596
597  DO il = 1, ncum
598    lwork(il) = .FALSE.
599  END DO
600
601! ---------------------------------------------------------------
602  DO i = minorig + 1, nl      !Loop on origin level "i"
603! ---------------------------------------------------------------
604
605    num1 = 0
606    DO il = 1, ncum
607      IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
608    END DO
609    IF (num1<=0) GO TO 789
610
611
612!JYG1    Find maximum of SIJ for J>I, if any.
613
614    Sx(:) = 0.
615
616    DO il = 1, ncum
617      IF (i>=icb(il) .AND. i<=inb(il)) THEN
618        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
619        Sbef(il) = max(0., signhpmh(il))
620      END IF
621    END DO
622
623    DO j = i + 1, nl
624      DO il = 1, ncum
625        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
626          IF (Sbef(il)<Sij(il,i,j)) THEN
627            Sx(il) = max(Sij(il,i,j), Sx(il))
628          END IF
629          Sbef(il) = Sij(il, i, j)
630        END IF
631      END DO
632    END DO
633
634
635    DO il = 1, ncum
636      IF (i>=icb(il) .AND. i<=inb(il)) THEN
637        lwork(il) = (nent(il,i)/=0)
638!!        rti = qnk(il) - ep(il, i)*clw(il, i)
639        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
640!jyg<
641        IF (cvflag_ice) THEN
642
643          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
644                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
645          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
646                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
647        ELSE
648
649          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
650                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
651          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
652                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
653        END IF
654!>jyg
655        IF (abs(denom)<0.01) denom = 0.01
656        Scrit(il) = min(anum/denom, 1.)
657        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
658
659!JYG1    Find new critical value Scrit2
660!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
661!                     Sij < Scrit2  => mixed draught will detrain at J>I
662
663        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
664                 Scrit(il)*max(0., signhpmh(il))
665
666        Scrit(il) = Scrit2
667
668!JYG    Correction pour la nouvelle logique; la correction pour ALT
669! est un peu au hazard
670        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
671        IF (alt<=0.0) Scrit(il) = 1.0
672
673        smax(il) = 0.0
674        ASij(il) = 0.0
675        sup(il) = 0.      ! upper S-value reached by descending draughts
676      END IF
677    END DO
678
679! ---------------------------------------------------------------
680    DO j = minorig, nl         !Loop on destination level "j"
681! ---------------------------------------------------------------
682
683      num2 = 0
684      DO il = 1, ncum
685        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
686            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
687            lwork(il)) num2 = num2 + 1
688      END DO
689      IF (num2<=0) GO TO 175
690
691! -----------------------------------------------
692      IF (j>i) THEN
693! -----------------------------------------------
694        DO il = 1, ncum
695          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
696              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
697              lwork(il)) THEN
698            IF (Sij(il,i,j)>0.0) THEN
699              Smid(il) = min(Sij(il,i,j), Scrit(il))
700              Sjmax(il) = Smid(il)
701              Sjmin(il) = Smid(il)
702              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
703                smin(il) = Smid(il)
704                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
705                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
706                Sjmin(il) = min(Sjmin(il), Scrit(il))
707                Sbef(il) = Sij(il, i, j)
708              END IF
709            END IF
710          END IF
711        END DO
712! -----------------------------------------------
713      ELSE IF (j==i) THEN
714! -----------------------------------------------
715        DO il = 1, ncum
716          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
717              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
718              lwork(il)) THEN
719            IF (Sij(il,i,j)>0.0) THEN
720              Smid(il) = 1.
721              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
722                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
723              Sjmin(il) = max(Sjmin(il), sup(il))
724              Sjmax(il) = 1.
725
726! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
727              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
728
729              smin(il) = 1.
730              Sbef(il) = max(0., signhpmh(il))
731              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
732            END IF
733          END IF
734        END DO
735! -----------------------------------------------
736      ELSE IF (j<i) THEN
737! -----------------------------------------------
738        DO il = 1, ncum
739          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
740              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
741              lwork(il)) THEN
742            IF (Sij(il,i,j)>0.0) THEN
743              Smid(il) = max(Sij(il,i,j), Scrit(il))
744              Sjmax(il) = Smid(il)
745              Sjmin(il) = Smid(il)
746              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
747                smax(il) = Smid(il)
748                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
749                Sjmax(il) = max(Sjmax(il), Scrit(il))
750                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
751                Sjmin(il) = max(Sjmin(il), Scrit(il))
752                Sbef(il) = Sij(il, i, j)
753              END IF
754              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
755                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
756            END IF
757          END IF
758        END DO
759! -----------------------------------------------
760      END IF
761! -----------------------------------------------
762
763
764      DO il = 1, ncum
765        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
766            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
767            lwork(il)) THEN
768          IF (Sij(il,i,j)>0.0) THEN
769!!            rti = qnk(il) - ep(il, i)*clw(il, i)
770            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
771            Qmixmax(il) = Qmix(Sjmax(il))
772            Qmixmin(il) = Qmix(Sjmin(il))
773            Rmixmax(il) = Rmix(Sjmax(il))
774            Rmixmin(il) = Rmix(Sjmin(il))
775            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
776            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
777
778            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
779
780! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
781            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
782              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
783            ELSE
784              Sigij(il, i, j) = 0.
785            END IF
786
787! --    Compute Qent, uent, vent according to the true mixing fraction
788            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
789            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
790            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)
791
792! --     Compute liquid water static energy of mixed draughts
793!    IF (j .GT. i) THEN
794!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
795!      awat=amax1(awat,0.0)
796!    ELSE
797!      awat = 0.
798!    ENDIF
799!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
800!    :         + Sigij(il,i,j)*H(il,i)
801!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
802!IM 301008 beg
803            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)
804
805!jyg<
806!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
807!            elij(il, i, j) = elij(il, i, j) + &
808!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
809!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
810!            elij(il, i, j) = elij(il, i, j) / &
811!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
812!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
813
814!       Computation of condensate amount Elij, taking into account the ice fraction frac
815!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
816!                 should be corrected.
817
818!  Heat capacity of mixed draught
819    cpm = cpd+Qent(il,i,j)*(cpv-cpd)
820
821    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
822            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
823            elij(il, i, j) = elij(il, i, j) + &
824                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
825                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
826            elij(il, i, j) = elij(il, i, j) / &
827                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
828                              (cpm*rrv*t(il,j)*t(il,j)))
829    ELSE
830            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
831            elij(il, i, j) = elij(il, i, j) + &
832                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
833                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
834            elij(il, i, j) = elij(il, i, j) / &
835                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
836                              (cpm*rrv*t(il,j)*t(il,j)))
837    ENDIF
838!>jyg
839            elij(il, i, j) = max(elij(il,i,j), 0.)
840
841            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
842
843            IF (j>i) THEN
844              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
845              awat = amax1(awat, 0.0)
846            ELSE
847              awat = 0.
848            END IF
849
850! PRINT *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
851! :         t(il,j))
852
853!jyg<
854!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
855! Mixed draught temperature at level j
856    IF (cvflag_ice .AND. frac(il,j) > 0.) THEN
857          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))
858          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
859    ELSE
860          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))
861          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
862    ENDIF
863!>jyg
864
865!IM 301008 end
866
867! PRINT *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
868! :               i,j,hent(il,i,j),Sigij(il,i,j)
869
870! --      ASij is the integral of P(F) over the relevant F interval
871            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
872                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
873
874          END IF
875        END IF
876      END DO ! il
877
878
879             
880#ifdef ISO
881      DO il=1,ncum
882          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
883          zfice(il) = MIN(MAX(zfice(il),0.0),1.0) 
884       IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
885           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
886           .AND. lwork(il) ) THEN
887        IF(sij(il,i,j).gt.0.0)THEN
888          DO ixt=1,ntraciso
889            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
890            xtent(ixt,il,i,j) = (1.-Sigij(il,i,j))*xtrti(ixt) &
891                    + Sigij(il,i,j)*xt(ixt,il,i)
892          enddo !do ixt = 1, ntraciso     
893#ifdef ISOVERIF
894#ifdef ISOTRAC         
895          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
896          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
897#endif         
898          IF (iso_eau.gt.0) THEN
899            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
900                qent(il,i,j),'cv3p_routines 718')   
901          endif
902#endif             
903         endif  !IF(sij(il,i,j).gt.0.0)THEN
904        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
905       enddo  !do il=1,ncum
906   
907#ifdef ISOVERIFIDRIS
908        DO il=1,ncum
909         CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 771')
910        enddo !do il=1,ncum   
911#endif     
912#ifdef ISOVERIF
913        IF (iso_eau.gt.0) THEN
914         DO il=1,ncum
915            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
916                qent(il,i,j),'cv3p_routines 744')   
917         enddo !do il=1,ncum   
918        endif   
919        DO il=1,ncum
920         CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
921                'cv3p_routines 749')
922        enddo !do il=1,ncum   
923        !WRITE(*,*) 'cv3p_mixing 788: CALL condiso_liq_ice_vectall'
924#endif     
925     
926       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
927                elij(1,i,j), &
928                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
929#ifdef ISOTRAC
930        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
931                elij(1,i,j), &
932                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
933#endif   
934        DO il=1,ncum
935        IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
936           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
937           .AND. lwork(il) ) THEN
938        IF(sij(il,i,j).gt.0.0)THEN
939         DO ixt = 1, ntraciso
940          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
941         enddo !do ixt = 1, ntraciso
942        endif !IF(sij(il,i,j).gt.0.0)THEN
943        endif !if ( i.ge.icb(il) .AND. i.le.inb(il) .AND.
944        enddo !do il=1,ncum
945
946#ifdef ISOVERIF 
947        IF (iso_eau.gt.0) THEN
948            DO il=1,ncum
949            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
950                qent(il,i,j),'cv3p_routines 744')
951            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
952                elij(il,i,j),'cv3p_routines 746')
953            enddo !do il=1,ncum   
954        endif
955#endif 
956!#ifdef ISOVERIF 
957
958#ifdef ISOTRAC   
959!        WRITE(*,*) 'cv3_routines tmp 1987,option_traceurs=',
960!     :           option_traceurs
961        IF (option_tmin.ge.1) THEN
962        DO il=1,ncum
963!        WRITE(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
964!     :           'tcond(il),rs(il,j)=',
965!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
966        ! colorier la vapeur résiduelle selon température de
967        ! condensation, et le condensat en un tag spécifique
968          IF ((elij(il,i,j).gt.0.0).AND.(qent(il,i,j).gt.0.0)) THEN
969            IF (option_traceurs.EQ.17) THEN
970             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
971                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
972                0.0,xtres, &
973                seuil_tag_tmin)
974            else !if (option_traceurs.EQ.17) THEN
975!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
976             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
977                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
978                seuil_tag_tmin)
979            endif !if (option_traceurs.EQ.17) THEN
980            DO ixt=1+niso,ntraciso
981               xtent(ixt,il,i,j)=xtres(ixt)
982            enddo     
983          endif !if (cond.gt.0.0) THEN
984        enddo !do il=1,ncum
985#ifdef ISOVERIF
986        DO il=1,ncum
987          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
988          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
989          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
990                'cv3_routines 2042')
991        enddo !do il=1,ncum
992#endif       
993        endif !if (option_tmin.ge.1) THEN
994#endif
995! END IF ISOTRAC
996
997#endif
998!!#ifdef ISO
999
1000!jyg!      DO k = 1, ntra
1001!jyg!        DO il = 1, ncum
1002!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1003!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1004!jyg!              lwork(il)) THEN
1005!jyg!            IF (Sij(il,i,j)>0.0) THEN
1006!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
1007!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
1008!jyg!            END IF
1009!jyg!          END IF
1010!jyg!        END DO
1011!jyg!      END DO
1012
1013! --    If I=J (detrainement and entrainement at the same level), then only the
1014! --    adiabatic ascent part of the mixture is considered
1015      IF (i==j) THEN
1016        DO il = 1, ncum
1017          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
1018              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
1019              lwork(il)) THEN
1020            IF (Sij(il,i,j)>0.0) THEN
1021!!              rti = qnk(il) - ep(il, i)*clw(il, i)
1022              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
1023!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
1024              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
1025                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
1026              Qent(il, i, i) = rti
1027              uent(il, i, i) = unk(il)
1028              vent(il, i, i) = vnk(il)
1029              hent(il, i, i) = hp(il, i)
1030              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1031              Sigij(il, i, i) = 0.
1032            END IF
1033          END IF
1034        END DO
1035
1036
1037#ifdef ISO
1038      IF (I .EQ. J) THEN
1039      DO il=1,ncum
1040          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
1041          zfice(il) = MIN(MAX(zfice(il),0.0),1.0) 
1042       IF ( i.ge.icb(il) .AND. i.le.inb(il) .AND. &
1043           j.ge.(icb(il)-1) .AND. j.le.inb(il) &
1044           .AND. lwork(il) ) THEN
1045        IF(sij(il,i,j).gt.0.0)THEN
1046          DO ixt=1,ntraciso
1047            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
1048            xtent(ixt,il,i,j) = xtrti(ixt)
1049          enddo !do ixt = 1, ntraciso
1050#ifdef ISOVERIFIDRIS
1051          CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 923')
1052#endif
1053#ifdef ISOVERIF 
1054#ifdef ISOTRAC         
1055          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
1056          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
1057#endif         
1058          IF (iso_eau.gt.0) THEN
1059            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
1060                qent(il,i,j),'cv3p_routines 718')   
1061          endif         
1062          CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
1063                'cv3p_routines 887')
1064#endif
1065         endif  !IF(sij(il,i,j).gt.0.0)THEN
1066        endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND.
1067       enddo  !do il=1,ncum
1068     
1069       !WRITE(*,*) 'cv3p_mixing 808: CALL condiso'
1070       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
1071                elij(1,i,j), &
1072                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1073
1074       DO il=1,ncum
1075         DO ixt = 1, ntraciso
1076          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
1077         enddo !do ixt = 1, ntraciso
1078#ifdef ISOVERIF
1079         IF (iso_eau.gt.0) THEN
1080           CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1081                elij(il,i,j),'cv3p_routines 1074')
1082         endif ! if (iso_eau.gt.0) THEN
1083#endif
1084        enddo !do il=1,ncum
1085       
1086#ifdef ISOTRAC
1087        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
1088                elij(1,i,j), &
1089                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
1090#ifdef ISOVERIF
1091        DO il=1,ncum
1092          CALL iso_verif_traceur(xt(1,il,i),'cv3p_routines 1967')
1093          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 1969')
1094        enddo !do il=1,ncum
1095#endif   
1096       
1097       IF (option_tmin.ge.1) THEN
1098        DO il=1,ncum
1099!        WRITE(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
1100!     :           'tcond(il),rs(il,j)=',
1101!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
1102        ! colorier la vapeur résiduelle selon température de
1103        ! condensation, et le condensat en un tag spécifique
1104          IF ((elij(il,i,j).gt.0.0).AND.(qent(il,i,j).gt.0.0)) THEN
1105            IF (option_traceurs.EQ.17) THEN
1106             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1107                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
1108                0.0,xtres, &
1109                seuil_tag_tmin)
1110            else !if (option_traceurs.EQ.17) THEN
1111!             WRITE(*,*) 'cv3 2002: il,i,j  =',il,i,j
1112             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
1113                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
1114                seuil_tag_tmin)
1115            endif !if (option_traceurs.EQ.17) THEN
1116            DO ixt=1+niso,ntraciso
1117               xtent(ixt,il,i,j)=xtres(ixt)
1118            enddo     
1119          endif !if (cond.gt.0.0) THEN
1120        enddo !do il=1,ncum
1121#ifdef ISOVERIF
1122        DO il=1,ncum
1123          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
1124          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
1125          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1126                'cv3_routines 2042')
1127        enddo !do il=1,ncum
1128#endif       
1129        endif !if (option_tmin.ge.1) THEN
1130!#ifdef ISOVERIF       
1131#endif
1132!#ifdef ISOTRAC       
1133#ifdef ISOVERIF       
1134        IF (iso_eau.gt.0) THEN
1135            DO il=1,ncum
1136            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
1137                qent(il,i,j),'cv3p_routines 883')
1138            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1139                elij(il,i,j),'cv3p_routines 885') 
1140            enddo !do il=1,ncum 
1141        endif !if (iso_eau.gt.0) THEN
1142#endif 
1143!#ifdef ISOVERIF           
1144        endif !IF (I .EQ. J) THEN
1145#endif 
1146!#ifdef ISO
1147
1148!jyg!        DO k = 1, ntra
1149!jyg!          DO il = 1, ncum
1150!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
1151!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
1152!jyg!                lwork(il)) THEN
1153!jyg!              IF (Sij(il,i,j)>0.0) THEN
1154!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
1155!jyg!              END IF
1156!jyg!            END IF
1157!jyg!          END DO
1158!jyg!        END DO
1159
1160      END IF
1161
1162! ---------------------------------------------------------------
1163175 END DO        ! End loop on destination level "j"
1164! ---------------------------------------------------------------
1165
1166    DO il = 1, ncum
1167      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
1168        ASij(il) = amax1(1.0E-16, ASij(il))
1169!jyg+lluis<
1170!!        ASij(il) = 1.0/ASij(il)
1171        ASij_inv(il) = 1.0/ASij(il)
1172!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
1173        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
1174!>jyg+lluis
1175        csum(il, i) = 0.0
1176      END IF
1177    END DO
1178
1179    DO j = minorig, nl
1180      DO il = 1, ncum
1181        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1182            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1183!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
1184          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
1185        END IF
1186      END DO
1187    END DO
1188
1189    DO j = minorig, nl
1190      DO il = 1, ncum
1191        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
1192            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
1193          csum(il, i) = csum(il, i) + Ment(il, i, j)
1194        END IF
1195      END DO
1196    END DO
1197
1198    DO il = 1, ncum
1199      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1200! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
1201        nent(il, i) = 0
1202! cc        Ment(il,i,i)=m(il,i)
1203        Ment(il, i, i) = 1.
1204!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
1205        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
1206        uent(il, i, i) = unk(il)
1207        vent(il, i, i) = vnk(il)
1208        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
1209
1210#ifdef ISO
1211      DO ixt = 1, ntraciso
1212        xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
1213        xtelij(ixt,il,i,i)=xtclw(ixt,il,i)*(1.-ep(il,i))
1214      enddo
1215#ifdef ISOVERIF
1216      IF (iso_eau.gt.0) THEN
1217        CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
1218                elij(il,i,j),'cv3p_routines 1204')
1219      endif ! if (iso_eau.gt.0) THEN
1220#endif
1221#endif
1222
1223#ifdef ISOTRAC         
1224        IF (option_tmin.ge.1) THEN
1225        ! colorier la vapeur résiduelle selon température de
1226        ! condensation, et le condensat en un tag spécifique
1227!        WRITE(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=',
1228!     :            il,i,j,xtent(:,il,i,j)
1229          IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN
1230            IF (option_traceurs.EQ.17) THEN
1231              CALL iso_recolorise_condensation(qent(il,i,i), &
1232                elij(il,i,i), &
1233                xt(1,il,1),xtclw(1,il,i),t(il,i),ep(il,i), &
1234                xtres, &
1235                seuil_tag_tmin)
1236            else !if (option_traceurs.EQ.17) THEN
1237              CALL iso_recolorise_condensation(qent(il,i,i), &
1238                elij(il,i,i), &
1239                xt(1,il,1),xtclw(1,il,i),rs(il,i),ep(il,i), &
1240                xtres, &
1241                seuil_tag_tmin)
1242            endif ! if (option_traceurs.EQ.17) THEN
1243            DO ixt=1+niso,ntraciso
1244              xtent(ixt,il,i,i)=xtres(ixt)
1245            enddo 
1246#ifdef ISOVERIF               
1247            DO ixt=1,niso
1248              CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
1249                'cv3p_mixing 2318',errmax,errmaxrel)
1250              CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
1251                'cv3p_mixing 2383')
1252            enddo
1253#endif               
1254          endif !if (cond.gt.0.0) THEN
1255#ifdef ISOVERIF         
1256          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
1257                qent(il,i,i),'cv3p_mixing 2321',errmax,errmaxrel)
1258          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), &
1259                elij(il,i,i),'cv3p_mixing 2321b',errmax,errmaxrel)
1260          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2322')
1261          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2323')
1262#endif       
1263        endif !if (option_tmin.ge.1) THEN
1264#endif
1265        IF (fl_cor_ebil .GE. 2) THEN
1266          hent(il, i, i) = hp(il,i)
1267          Sigij(il, i, i) = 0.0
1268        ELSE
1269          Sij(il, i, i) = 0.0
1270        ENDIF
1271      END IF
1272    END DO ! il
1273
1274!jyg!    DO j = 1, ntra
1275!jyg!      DO il = 1, ncum
1276!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
1277!jyg!! cc     :     .AND. csum(il,i).lt.m(il,i) ) THEN
1278!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
1279!jyg!        END IF
1280!jyg!      END DO
1281!jyg!    END DO
1282
1283! ---------------------------------------------------------------
1284789 END DO              ! End loop on origin level "i"
1285! ---------------------------------------------------------------
1286
1287#ifdef ISO
1288#ifdef ISOVERIF
1289!       WRITE(*,*) 'cv3p_mixing 2540: ', &
1290!        'verif finale en sortant de cv3p_mixing'
1291!       WRITE(*,*) 'qent,xtent(1,1,1)=',qent(1,1,1),xtent(iso_eau,1,1,1)
1292       DO im = 1, nd
1293       DO jm = 1, nd
1294        DO il = 1, ncum
1295          IF (iso_eau.gt.0) THEN
1296            CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1297              elij(il,im,jm),'cv3p_mixing 2110',errmax,errmaxrel)
1298            CALL iso_verif_egalite_choix(xtent(iso_eau,il,im,jm), &
1299              qent(il,im,jm),'cv3p_mixing 2112',errmax,errmaxrel)
1300          endif !if (iso_eau>0) THEN
1301#ifdef ISOTRAC
1302        CALL iso_verif_traceur_justmass(xtelij(1,il,im,jm), &
1303                       'cv3p_routine 2250')
1304#endif           
1305        enddo !do il = 1, nloc
1306       enddo !do jm = 1, klev
1307       enddo !do im = 1, klev
1308#endif
1309#endif 
1310
1311#ifdef ISO
1312#ifdef ISOTRAC
1313        ! seulement à la fin on taggue le condensat
1314        IF (option_cond.ge.1) THEN
1315         DO im = 1, nd
1316         DO jm = 1, nd
1317         DO il = 1, ncum
1318           ! colorier le condensat en un tag spécifique
1319           DO ixt=niso+1,ntraciso
1320             IF (index_zone(ixt).EQ.izone_cond) THEN
1321                xtelij(ixt,il,im,jm)=xtelij(index_iso(ixt),il,im,jm)
1322             else !if (index_zone(ixt).EQ.izone_cond) THEN
1323                xtelij(ixt,il,im,jm)=0.0
1324             endif !if (index_zone(ixt).EQ.izone_cond) THEN
1325           enddo !do ixt=1,ntraciso     
1326#ifdef ISOVERIF
1327        CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
1328                elij(il,im,jm),'cv3p_mixing 2408',errmax,errmaxrel)
1329        CALL iso_verif_traceur(xtelij(1,il,im,jm), &
1330               'cv3p_mixing 358')
1331#endif     
1332         enddo !do il = 1, ncum   
1333         enddo !do jm = 1, nd
1334         enddo !do im = 1, nd
1335        ! 17 juin 2020: on supprime ci-dessous car normalement xtclw a déjà été
1336        ! recolorié ds cv3_undilute
1337!         do im = 1, nd
1338!         do il = 1, ncum   
1339!           ! colorier le condensat en un tag spécifique
1340!           do ixt=niso+1,ntraciso
1341!             if (index_zone(ixt).EQ.izone_cond) THEN
1342!                xtclw(ixt,il,im)=xtclw(index_iso(ixt),il,im)
1343!             else !if (index_zone(ixt).EQ.izone_cond) THEN
1344!                xtclw(ixt,il,im)=0.0
1345!             endif !if (index_zone(ixt).EQ.izone_cond) THEN
1346!           enddo !do ixt=1,ntraciso     
1347!#ifdef ISOVERIF
1348!        CALL iso_verif_egalite_choix(xtclw(iso_eau,il,im), &
1349!                clw(il,im),'cv3p_mixing 2427',errmax,errmaxrel)
1350!        CALL iso_verif_traceur(xtclw(1,il,im), &
1351!               'cv3p_mixing 358')
1352!        if (iso_verif_positif_nostop(xtclw(itZonIso( &
1353!                izone_cond,iso_eau),i,k)-xtclw(iso_eau,i,k) &
1354!                ,'cv3p_mixing 909').EQ.1) THEN
1355!               WRITE(*,*) 'i,k=',i,k
1356!               WRITE(*,*) 'xtclw=',xtclw(:,i,k)
1357!               WRITE(*,*) 'niso,ntraciso,index_zone,izone_cond=', &
1358!                  niso,ntraciso,index_zone,izone_cond     
1359!               stop
1360!         endif !if (iso_verif_positif_nostop(xtclw(itZonIso(
1361!#endif             
1362!         enddo !do il = 1, ncum   
1363!         enddo !do im = 1, nd
1364!         WRITE(*,*) 'xtclw(:,1,2)=',xtclw(:,1,2)
1365        endif !if (option_tmin.EQ.1) THEN
1366#endif
1367#endif         
1368
1369
1370
1371END SUBROUTINE cv3p_mixing
1372
Note: See TracBrowser for help on using the repository browser.