source: LMDZ4/trunk/libf/phylmd/cv3p_mixing.F @ 880

Last change on this file since 880 was 879, checked in by Laurent Fairhead, 17 years ago

Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ...
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.4 KB
Line 
1      SUBROUTINE cv3p_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
2     :                    ,ph,t,rr,rs,u,v,tra,h,lv,qnk
3     :                    ,unk,vnk,hp,tv,tvp,ep,clw,sig
4     :                    ,ment,qent,hent,uent,vent
5     :                   ,sij,elij,supmax,ments,qents,traent)
6***************************************************************
7*                                                             *
8* CV3P_MIXING : compute mixed draught properties and,         *
9*               within a scaling factor, mixed draught        *
10*               mass fluxes.                                  *
11* written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
12* modified by :                                               *
13***************************************************************
14*
15      implicit none
16c
17#include "cvthermo.h"
18#include "cv3param.h"
19#include "YOMCST2.h"
20
21c inputs:
22      integer ncum, nd, na, ntra, nloc
23      integer icb(nloc), inb(nloc), nk(nloc)
24      real sig(nloc,nd)
25      real qnk(nloc),unk(nloc),vnk(nloc)
26      real ph(nloc,nd+1)
27      real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
28      real u(nloc,nd), v(nloc,nd)
29      real tra(nloc,nd,ntra) ! input of convect3
30      real lv(nloc,na)
31      real h(nloc,na)  !liquid water static energy of environment
32      real hp(nloc,na) !liquid water static energy of air shed from adiab. asc.
33      real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
34
35c outputs:
36      real ment(nloc,na,na), qent(nloc,na,na)
37      real uent(nloc,na,na), vent(nloc,na,na)
38      real sij(nloc,na,na), elij(nloc,na,na)
39      real supmax(nloc,na)     ! Highest mixing fraction of mixed updraughts
40                               ! with the sign of (h-hp)
41      real traent(nloc,nd,nd,ntra)
42      real ments(nloc,nd,nd), qents(nloc,nd,nd)
43      real sigij(nloc,nd,nd)
44      real hent(nloc,nd,nd)
45
46c local variables:
47      integer i, j, k, il, im, jm
48      integer num1, num2
49      integer nent(nloc,na)
50      real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
51      real alt, delp, delm
52      real Qmixmax(nloc), Rmixmax(nloc), SQmRmax(nloc)
53      real Qmixmin(nloc), Rmixmin(nloc), SQmRmin(nloc)
54      real signhpmh(nloc)
55      real Sx, Scrit2
56      integer Jx
57      real smid(nloc), sjmin(nloc), sjmax(nloc)
58      real Sbef(nloc), Sup(nloc), Smin(nloc)
59      real asij(nloc), smax(nloc), scrit(nloc)
60      real csum(nloc,nd)
61      real awat
62      logical lwork(nloc)
63c
64      REAL amxupcrit, df, ff
65      INTEGER nstep
66C
67c --   Mixing probability distribution functions
68c
69      real Qcoef1,Qcoef2,QFF,QFFF,Qmix,Rmix,Qmix1,Rmix1,Qmix2,Rmix2,F
70      Qcoef1(F) = tanh(F/gammas)
71      Qcoef2(F) = ( tanh(F/gammas) + gammas *
72     $            log(cosh((1.- F)/gammas)/cosh(F/gammas)))
73      QFF(F) = Max(Min(F,1.),0.)
74      QFFF(F) = Min(QFF(F),scut)
75      Qmix1(F) = ( tanh((QFF(F) - Fmax)/gammas)+Qcoef1max )/
76     $           Qcoef2max
77      Rmix1(F) = ( gammas*log(cosh((QFF(F)-Fmax)/gammas))
78     1             +QFF(F)*Qcoef1max ) / Qcoef2max
79      Qmix2(F) = -Log(1.-QFFF(F))/scut
80      Rmix2(F) = (QFFF(F)+(1.-QFF(F))*Log(1.-QFFF(F)))/scut
81      Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
82      Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)
83C
84      INTEGER ifrst
85      DATA ifrst/0/
86C
87
88c=====================================================================
89c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
90c=====================================================================
91c
92c -- Initialize mixing PDF coefficients
93      IF (ifrst .EQ. 0) THEN
94        ifrst = 1
95        Qcoef1max = Qcoef1(Fmax)
96        Qcoef2max = Qcoef2(Fmax)
97c
98      ENDIF
99c
100
101c ori        do 360 i=1,ncum*nlp
102        do 361 j=1,nl
103        do 360 i=1,ncum
104          nent(i,j)=0
105c in convect3, m is computed in cv3_closure
106c ori          m(i,1)=0.0
107 360    continue
108 361    continue
109
110c ori      do 400 k=1,nlp
111c ori       do 390 j=1,nlp
112      do 400 j=1,nl
113       do 390 k=1,nl
114          do 385 i=1,ncum
115            qent(i,k,j)=rr(i,j)
116            uent(i,k,j)=u(i,j)
117            vent(i,k,j)=v(i,j)
118            elij(i,k,j)=0.0
119            hent(i,k,j)=0.0
120            ment(i,k,j)=0.0
121            sij(i,k,j)=0.0
122 385      continue
123 390    continue
124 400  continue
125
126      do k=1,ntra
127       do j=1,nd  ! instead nlp
128        do i=1,nd ! instead nlp
129         do il=1,ncum
130            traent(il,i,j,k)=tra(il,j,k)
131         enddo
132        enddo
133       enddo
134      enddo
135
136c=====================================================================
137c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
138c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
139c --- FRACTION (sij)
140c=====================================================================
141
142      do 750 i=minorig+1, nl
143
144       do 710 j=minorig,nl
145        do 700 il=1,ncum
146         if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
147     :      (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
148
149          rti=qnk(il)-ep(il,i)*clw(il,i)
150          bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
151          anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
152          denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
153          dei=denom
154          if(abs(dei).lt.0.01)dei=0.01
155          sij(il,i,j)=anum/dei
156          sij(il,i,i)=1.0
157          altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
158          altem=altem/bf2
159          cwat=clw(il,j)*(1.-ep(il,j))
160          stemp=sij(il,i,j)
161          if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
162     :                 .and.j.gt.i)then
163           anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
164           denom=denom+lv(il,j)*(rr(il,i)-rti)
165           if(abs(denom).lt.0.01)denom=0.01
166           sij(il,i,j)=anum/denom
167           altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
168           altem=altem-(bf2-1.)*cwat
169          end if
170              if(sij(il,i,j).gt.0.0)then
171ccc                 ment(il,i,j)=m(il,i)
172                 ment(il,i,j)=1.
173                 elij(il,i,j)=altem
174                 elij(il,i,j)=amax1(0.0,elij(il,i,j))
175                 nent(il,i)=nent(il,i)+1
176              endif
177
178         sij(il,i,j)=amax1(0.0,sij(il,i,j))
179         sij(il,i,j)=amin1(1.0,sij(il,i,j))
180         endif ! new
181 700   continue
182 710  continue
183
184c
185c   ***   if no air can entrain at level i assume that updraft detrains  ***
186c   ***   at that level and calculate detrained air flux and properties  ***
187c
188
189c@      do 170 i=icb(il),inb(il)
190
191      do 740 il=1,ncum
192      if ((i.ge.icb(il)).and.(i.le.inb(il))
193     :                  .and.(nent(il,i).eq.0)) then
194c@      if(nent(il,i).eq.0)then
195ccc      ment(il,i,i)=m(il,i)
196      ment(il,i,i)=1.
197      qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
198      uent(il,i,i)=unk(il)
199      vent(il,i,i)=vnk(il)
200      elij(il,i,i)=clw(il,i)*(1.-ep(il,i))
201      sij(il,i,i)=0.0
202      end if
203 740  continue
204 750  continue
205
206      do j=1,ntra
207       do i=minorig+1,nl
208        do il=1,ncum
209         if (i.ge.icb(il) .and. i.le.inb(il)
210     :                    .and. nent(il,i).eq.0) then
211          traent(il,i,i,j)=tra(il,nk(il),j)
212         endif
213        enddo
214       enddo
215      enddo
216
217      do 100 j=minorig,nl
218      do 101 i=minorig,nl
219      do 102 il=1,ncum
220      if ((j.ge.(icb(il)-1)).and.(j.le.inb(il))
221     :    .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
222       sigij(il,i,j)=sij(il,i,j)
223      endif
224 102  continue
225 101  continue
226 100  continue
227c@      enddo
228
229c@170   continue
230
231c=====================================================================
232c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
233c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
234c=====================================================================
235
236      call zilch(csum,nloc*nd)
237
238      do il=1,ncum
239       lwork(il) = .FALSE.
240      enddo
241
242      DO 789 i=minorig+1,nl
243
244      num1=0
245      do il=1,ncum
246       if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
247      enddo
248      if (num1.le.0) goto 789
249
250
251      do 781 il=1,ncum
252       if ( i.ge.icb(il) .and. i.le.inb(il) ) then
253        lwork(il)=(nent(il,i).ne.0)
254        Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
255        qp=qnk(il)-ep(il,i)*clw(il,i)
256        anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
257     :           +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
258        denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
259     :           +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
260        if(abs(denom).lt.0.01)denom=0.01
261        scrit(il)=min(anum/denom,1.)
262        alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
263c
264cjyg1    Find maximum of SIJ for J>I, if any, and new critical value Scrit2
265c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
266c                    Sij < Scrit2  => mixed draught will detrain at J>I
267c
268       Sx = 0.
269       Jx = 0.
270       Sbef(il) = max(0.,signhpmh(il))
271       DO j = i+1,inb(il)
272         IF (Sbef(il) .LT. Sij(il,i,j)) THEN
273           Sx = max(Sij(il,i,j),Sx)
274           Jx = J
275         ENDIF
276         Sbef(il) = Sij(il,i,j)
277       ENDDO
278c
279       Scrit2 = min(Scrit(il),Sx)*max(0.,-signhpmh(il))
280     :         +Scrit(il)*max(0.,signhpmh(il))
281c
282       Scrit(il) = Scrit2
283c
284cjyg    Correction pour la nouvelle logique; la correction pour ALT
285c       est un peu au hazard
286       if(scrit(il).le.0.0)scrit(il)=0.0
287       if(alt.le.0.0) scrit(il)=1.0
288C
289        smax(il)=0.0
290        asij(il)=0.0
291       Sup(il)=0.     ! upper S-value reached by descending draughts
292       endif
293781   continue
294
295      do 175 j=minorig,nl
296
297      num2=0
298      do il=1,ncum
299       if ( i.ge.icb(il) .and. i.le.inb(il) .and.
300     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
301     :      .and. lwork(il) ) num2=num2+1
302      enddo
303      if (num2.le.0) goto 175
304
305c -----------------------------------------------
306         IF (j .GT. i) THEN
307c -----------------------------------------------
308      do 782 il=1,ncum
309      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
310     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
311     :      .and. lwork(il) ) then
312       if(sij(il,i,j).gt.0.0)then
313           Smid(il)=min(Sij(il,i,j),Scrit(il))
314           Sjmax(il)=Smid(il)
315           Sjmin(il)=Smid(il)
316           IF (Smid(il) .LT. Smin(il) .AND.
317     1                         Sij(il,i,j+1) .LT. Smid(il)) THEN
318             Smin(il)=Smid(il)
319             Sjmax(il)=min( (Sij(il,i,j+1)+Sij(il,i,j))/2. ,
320     1                  Sij(il,i,j) ,
321     1                  Scrit(il) )
322             Sjmin(il)=max( (Sbef(il)+Sij(il,i,j))/2. ,
323     1                  Sij(il,i,j) )
324             Sjmin(il)=min(Sjmin(il),Scrit(il))
325             Sbef(il) = Sij(il,i,j)
326           ENDIF
327      endif
328      endif
329782   continue
330c -----------------------------------------------
331         ELSE IF (j .EQ. i) THEN
332c -----------------------------------------------
333      do 783 il=1,ncum
334      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
335     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
336     :      .and. lwork(il) ) then
337       if(sij(il,i,j).gt.0.0)then
338           Smid(il) = 1.
339           Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2.,Scrit(il))
340     1                                         *max(0.,-signhpmh(il))
341     1            +min((Sij(il,i,j+1)+Smid(il))/2.,Scrit(il))
342     1                                         *max(0., signhpmh(il))
343           Sjmin(il) = max(Sjmin(il),Sup(il))
344           Sjmax(il) = 1.
345c
346c-           preparation des variables Scrit, Smin et Sbef pour la partie j>i
347           Scrit(il) = min(Sjmin(il),Sjmax(il),Scrit(il))
348
349           Smin(il) = 1.
350           Sbef(il) = max(0.,signhpmh(il))
351           Supmax(il,i) = sign(Scrit(il),-signhpmh(il))
352      endif
353      endif
354783   continue
355c -----------------------------------------------
356         ELSE IF ( j .LT. i) THEN
357c -----------------------------------------------
358      do 784 il=1,ncum
359      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
360     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
361     :      .and. lwork(il) ) then
362       if(sij(il,i,j).gt.0.0)then
363           Smid(il)=max(Sij(il,i,j),Scrit(il))
364           Sjmax(il) = Smid(il)
365           Sjmin(il) = Smid(il)
366           IF (Smid(il) .GT. Smax(il) .AND.
367     1                          Sij(il,i,j+1) .GT. Smid(il)) THEN
368             Smax(il) = Smid(il)
369             Sjmax(il) = max( (Sij(il,i,j+1)+Sij(il,i,j))/2. ,
370     1                                               Sij(il,i,j) )
371             Sjmax(il) = max(Sjmax(il),Scrit(il))
372             Sjmin(il) = min( (Sbef(il)+Sij(il,i,j))/2. ,
373     1                                               Sij(il,i,j) )
374             Sjmin(il) = max(Sjmin(il),Scrit(il))
375             Sbef(il) = Sij(il,i,j)
376           ENDIF
377          IF (abs(Sjmin(il)-Sjmax(il)) .GT. 1.e-10) Sup(il)=
378     1                            max(Sjmin(il),Sjmax(il),Sup(il))
379      endif
380      endif
381784   continue
382c -----------------------------------------------
383         END IF
384c -----------------------------------------------
385c
386c
387      do il=1,ncum
388      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
389     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
390     :      .and. lwork(il) ) then
391       if(sij(il,i,j).gt.0.0)then
392          rti=rr(il,1)-ep(il,i)*clw(il,i)
393         Qmixmax(il)=Qmix(Sjmax(il))
394         Qmixmin(il)=Qmix(Sjmin(il))
395         Rmixmax(il)=Rmix(Sjmax(il))
396         Rmixmin(il)=Rmix(Sjmin(il))
397         SQmRmax(il)= Sjmax(il)*Qmix(Sjmax(il))-Rmix(Sjmax(il))
398         SQmRmin(il)= Sjmin(il)*Qmix(Sjmin(il))-Rmix(Sjmin(il))
399c
400         Ment(il,i,j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il,i,j)
401c
402c    Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
403         IF (abs(Qmixmax(il)-Qmixmin(il)) .GT. 1.e-10) THEN
404           Sigij(il,i,j) =
405     :           (SQmRmax(il)-SQmRmin(il))/(Qmixmax(il)-Qmixmin(il))
406         ELSE
407           Sigij(il,i,j) = 0.
408         ENDIF
409c
410c --    Compute Qent, uent, vent according to the true mixing fraction
411        Qent(il,i,j) = (1.-Sigij(il,i,j))*rti
412     :               + Sigij(il,i,j)*rr(il,i)
413        uent(il,i,j) = (1.-Sigij(il,i,j))*unk(il)
414     :               + Sigij(il,i,j)*u(il,i)
415        vent(il,i,j) = (1.-Sigij(il,i,j))*vnk(il)
416     :               + Sigij(il,i,j)*v(il,i)
417c
418c--     Compute liquid water static energy of mixed draughts
419        IF (j .GT. i) THEN
420         awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
421         awat=amax1(awat,0.0)
422        ELSE
423         awat = 0.
424        ENDIF
425        Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
426     :         + Sigij(il,i,j)*H(il,i)
427     :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
428c
429c      print *,'mix : i,j,hent(il,i,j),sigij(il,i,j) ',
430c     :               i,j,hent(il,i,j),sigij(il,i,j)
431c
432c --      ASij is the integral of P(F) over the relevant F interval
433         ASij(il) = ASij(il)
434     1               + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il)
435     1                    -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
436c
437      endif
438      endif
439      enddo
440       do k=1,ntra
441         do il=1,ncum
442          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
443     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il))
444     :      .and. lwork(il) ) then
445          if(sij(il,i,j).gt.0.0)then
446            traent(il,i,j,k)=sigij(il,i,j)*tra(il,i,k)
447     :            +(1.-sigij(il,i,j))*tra(il,nk(il),k)
448          endif
449          endif
450         enddo
451       enddo
452c
453c --    If I=J (detrainement and entrainement at the same level), then only the
454c --    adiabatic ascent part of the mixture is considered
455        IF (I .EQ. J) THEN
456      do il=1,ncum
457      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
458     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
459     :      .and. lwork(il) ) then
460       if(sij(il,i,j).gt.0.0)then
461          rti=qnk(il)-ep(il,i)*clw(il,i)
462ccc          Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
463          Ment(il,i,i) = abs(Qmixmax(il)*(1.-Sjmax(il))
464     1                    +Rmixmax(il)
465     1                    -Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
466          Qent(il,i,i) = rti
467          uent(il,i,i) = unk(il)
468          vent(il,i,i) = vnk(il)
469          Hent(il,i,i) = hp(il,i)
470          Elij(il,i,i) = clw(il,i)*(1.-ep(il,i))
471          Sigij(il,i,i) = 0.
472      endif
473      endif
474      enddo
475       do k=1,ntra
476         do il=1,ncum
477          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
478     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il))
479     :      .and. lwork(il) ) then
480          if(sij(il,i,j).gt.0.0)then
481            traent(il,i,i,k)=tra(il,nk(il),k)
482          endif
483          endif
484         enddo
485       enddo
486c
487        ENDIF
488c
489175   continue
490
491      do il=1,ncum
492       if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
493        asij(il)=amax1(1.0e-16,asij(il))
494        asij(il)=1.0/asij(il)
495        csum(il,i)=0.0
496       endif
497      enddo
498
499      do 180 j=minorig,nl
500       do il=1,ncum
501        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
502     :   .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
503         ment(il,i,j)=ment(il,i,j)*asij(il)
504        endif
505       enddo
506180   continue
507
508      do 197 j=minorig,nl
509       do il=1,ncum
510        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
511     :   .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
512         csum(il,i)=csum(il,i)+ment(il,i,j)
513        endif
514       enddo
515197   continue
516
517      do il=1,ncum
518       if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
519     :     .and. csum(il,i).lt.1. ) then
520ccc     :     .and. csum(il,i).lt.m(il,i) ) then
521        nent(il,i)=0
522ccc        ment(il,i,i)=m(il,i)
523        ment(il,i,i)=1.
524        qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
525        uent(il,i,i)=unk(il)
526        vent(il,i,i)=vnk(il)
527        elij(il,i,i)=clw(il,i)*(1.-ep(il,i))
528        sij(il,i,i)=0.0
529       endif
530      enddo ! il
531
532      do j=1,ntra
533       do il=1,ncum
534        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
535     :     .and. csum(il,i).lt.1. ) then
536ccc     :     .and. csum(il,i).lt.m(il,i) ) then
537         traent(il,i,i,j)=tra(il,nk(il),j)
538        endif
539       enddo
540      enddo
541c
542789   continue
543c
544      return
545      end
546
Note: See TracBrowser for help on using the repository browser.