/[lmdze]/trunk/Sources/phylmd/CV30_routines/cv30_mixing.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/CV30_routines/cv30_mixing.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/phylmd/CV3_routines/cv3_mixing.f
File size: 10514 byte(s)
Changed all ".f90" suffixes to ".f".
1 guez 47
2     SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb &
3     ,ph,t,rr,rs,u,v,tra,h,lv,qnk &
4     ,hp,tv,tvp,ep,clw,m,sig &
5     ,ment,qent,uent,vent, nent, sij,elij,ments,qents,traent)
6 guez 69 use cv3_param_m
7 guez 47 use cvthermo
8     implicit none
9    
10     !---------------------------------------------------------------------
11     ! a faire:
12     ! - changer rr(il,1) -> qnk(il)
13     ! - vectorisation de la partie normalisation des flux (do 789...)
14     !---------------------------------------------------------------------
15    
16    
17     ! inputs:
18 guez 69 integer, intent(in):: ncum, nd, na, ntra, nloc
19 guez 47 integer icb(nloc), inb(nloc), nk(nloc)
20     real sig(nloc,nd)
21     real qnk(nloc)
22     real ph(nloc,nd+1)
23     real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
24     real u(nloc,nd), v(nloc,nd)
25     real tra(nloc,nd,ntra) ! input of convect3
26     real lv(nloc,na), h(nloc,na), hp(nloc,na)
27     real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
28     real m(nloc,na) ! input of convect3
29    
30     ! outputs:
31     real ment(nloc,na,na), qent(nloc,na,na)
32     real uent(nloc,na,na), vent(nloc,na,na)
33     real sij(nloc,na,na), elij(nloc,na,na)
34     real traent(nloc,nd,nd,ntra)
35     real ments(nloc,nd,nd), qents(nloc,nd,nd)
36     real sigij(nloc,nd,nd)
37     integer nent(nloc,nd)
38    
39     ! local variables:
40     integer i, j, k, il, im, jm
41     integer num1, num2
42     real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
43     real alt, smid, sjmin, sjmax, delp, delm
44     real asij(nloc), smax(nloc), scrit(nloc)
45     real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
46     real wgh
47     real zm(nloc,na)
48     logical lwork(nloc)
49    
50     !=====================================================================
51     ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
52     !=====================================================================
53    
54     do 361 j=1,nl
55     do 360 i=1,ncum
56     nent(i,j)=0
57     ! in convect3, m is computed in cv3_closure
58     ! ori m(i,1)=0.0
59     360 continue
60     361 continue
61    
62     ! ori do 400 k=1,nlp
63     ! ori do 390 j=1,nlp
64     do 400 j=1,nl
65     do 390 k=1,nl
66     do 385 i=1,ncum
67     qent(i,k,j)=rr(i,j)
68     uent(i,k,j)=u(i,j)
69     vent(i,k,j)=v(i,j)
70     elij(i,k,j)=0.0
71     !ym ment(i,k,j)=0.0
72     !ym sij(i,k,j)=0.0
73     385 continue
74     390 continue
75     400 continue
76    
77     !ym
78     ment(1:ncum,1:nd,1:nd)=0.0
79     sij(1:ncum,1:nd,1:nd)=0.0
80    
81     zm(:,:)=0.
82    
83     !=====================================================================
84     ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
85     ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
86     ! --- FRACTION (sij)
87     !=====================================================================
88    
89     do 750 i=minorig+1, nl
90    
91     do 710 j=minorig,nl
92     do 700 il=1,ncum
93     if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
94     (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
95    
96     rti=rr(il,1)-ep(il,i)*clw(il,i)
97     bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
98     anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
99     denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
100     dei=denom
101     if(abs(dei).lt.0.01)dei=0.01
102     sij(il,i,j)=anum/dei
103     sij(il,i,i)=1.0
104     altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
105     altem=altem/bf2
106     cwat=clw(il,j)*(1.-ep(il,j))
107     stemp=sij(il,i,j)
108     if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &
109     .and.j.gt.i)then
110     anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
111     denom=denom+lv(il,j)*(rr(il,i)-rti)
112     if(abs(denom).lt.0.01)denom=0.01
113     sij(il,i,j)=anum/denom
114     altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
115     altem=altem-(bf2-1.)*cwat
116     end if
117     if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
118     qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti
119     uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))
120     vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))
121     !!!! do k=1,ntra
122     !!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
123     !!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
124     !!!! end do
125     elij(il,i,j)=altem
126     elij(il,i,j)=amax1(0.0,elij(il,i,j))
127     ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
128     nent(il,i)=nent(il,i)+1
129     end if
130     sij(il,i,j)=amax1(0.0,sij(il,i,j))
131     sij(il,i,j)=amin1(1.0,sij(il,i,j))
132     endif ! new
133     700 continue
134     710 continue
135    
136     !
137     ! *** if no air can entrain at level i assume that updraft detrains ***
138     ! *** at that level and calculate detrained air flux and properties ***
139     !
140    
141     !@ do 170 i=icb(il),inb(il)
142    
143     do 740 il=1,ncum
144     if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
145     !@ if(nent(il,i).eq.0)then
146     ment(il,i,i)=m(il,i)
147     qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
148     uent(il,i,i)=u(il,nk(il))
149     vent(il,i,i)=v(il,nk(il))
150     elij(il,i,i)=clw(il,i)
151     !MAF sij(il,i,i)=1.0
152     sij(il,i,i)=0.0
153     end if
154     740 continue
155     750 continue
156    
157     do 100 j=minorig,nl
158     do 101 i=minorig,nl
159     do 102 il=1,ncum
160     if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
161     .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
162     sigij(il,i,j)=sij(il,i,j)
163     endif
164     102 continue
165     101 continue
166     100 continue
167     !@ enddo
168    
169     !@170 continue
170    
171     !=====================================================================
172     ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
173     ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
174     !=====================================================================
175    
176     call zilch(asum,nloc*nd)
177     call zilch(csum,nloc*nd)
178     call zilch(csum,nloc*nd)
179    
180     do il=1,ncum
181     lwork(il) = .FALSE.
182     enddo
183    
184     DO 789 i=minorig+1,nl
185    
186     num1=0
187     do il=1,ncum
188     if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
189     enddo
190     if (num1.le.0) goto 789
191    
192    
193     do 781 il=1,ncum
194     if ( i.ge.icb(il) .and. i.le.inb(il) ) then
195     lwork(il)=(nent(il,i).ne.0)
196     qp=rr(il,1)-ep(il,i)*clw(il,i)
197     anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) &
198     +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
199     denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) &
200     +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
201     if(abs(denom).lt.0.01)denom=0.01
202     scrit(il)=anum/denom
203     alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
204     if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
205     smax(il)=0.0
206     asij(il)=0.0
207     endif
208     781 continue
209    
210     do 175 j=nl,minorig,-1
211    
212     num2=0
213     do il=1,ncum
214     if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
215     j.ge.(icb(il)-1) .and. j.le.inb(il) &
216     .and. lwork(il) ) num2=num2+1
217     enddo
218     if (num2.le.0) goto 175
219    
220     do 782 il=1,ncum
221     if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
222     j.ge.(icb(il)-1) .and. j.le.inb(il) &
223     .and. lwork(il) ) then
224    
225     if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
226     wgh=1.0
227     if(j.gt.i)then
228     sjmax=amax1(sij(il,i,j+1),smax(il))
229     sjmax=amin1(sjmax,scrit(il))
230     smax(il)=amax1(sij(il,i,j),smax(il))
231     sjmin=amax1(sij(il,i,j-1),smax(il))
232     sjmin=amin1(sjmin,scrit(il))
233     if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
234     smid=amin1(sij(il,i,j),scrit(il))
235     else
236     sjmax=amax1(sij(il,i,j+1),scrit(il))
237     smid=amax1(sij(il,i,j),scrit(il))
238     sjmin=0.0
239     if(j.gt.1)sjmin=sij(il,i,j-1)
240     sjmin=amax1(sjmin,scrit(il))
241     endif
242     delp=abs(sjmax-smid)
243     delm=abs(sjmin-smid)
244     asij(il)=asij(il)+wgh*(delp+delm)
245     ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh
246     endif
247     endif
248     782 continue
249    
250     175 continue
251    
252     do il=1,ncum
253     if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
254     asij(il)=amax1(1.0e-16,asij(il))
255     asij(il)=1.0/asij(il)
256     asum(il,i)=0.0
257     bsum(il,i)=0.0
258     csum(il,i)=0.0
259     endif
260     enddo
261    
262     do 180 j=minorig,nl
263     do il=1,ncum
264     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
265     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
266     ment(il,i,j)=ment(il,i,j)*asij(il)
267     endif
268     enddo
269     180 continue
270    
271     do 190 j=minorig,nl
272     do il=1,ncum
273     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
274     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
275     asum(il,i)=asum(il,i)+ment(il,i,j)
276     ment(il,i,j)=ment(il,i,j)*sig(il,j)
277     bsum(il,i)=bsum(il,i)+ment(il,i,j)
278     endif
279     enddo
280     190 continue
281    
282     do il=1,ncum
283     if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
284     bsum(il,i)=amax1(bsum(il,i),1.0e-16)
285     bsum(il,i)=1.0/bsum(il,i)
286     endif
287     enddo
288    
289     do 195 j=minorig,nl
290     do il=1,ncum
291     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
292     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
293     ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i)
294     endif
295     enddo
296     195 continue
297    
298     do 197 j=minorig,nl
299     do il=1,ncum
300     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
301     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
302     csum(il,i)=csum(il,i)+ment(il,i,j)
303     endif
304     enddo
305     197 continue
306    
307     do il=1,ncum
308     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
309     .and. csum(il,i).lt.m(il,i) ) then
310     nent(il,i)=0
311     ment(il,i,i)=m(il,i)
312     qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
313     uent(il,i,i)=u(il,nk(il))
314     vent(il,i,i)=v(il,nk(il))
315     elij(il,i,i)=clw(il,i)
316     !MAF sij(il,i,i)=1.0
317     sij(il,i,i)=0.0
318     endif
319     enddo ! il
320    
321     789 continue
322     !
323     ! MAF: renormalisation de MENT
324     do jm=1,nd
325     do im=1,nd
326     do il=1,ncum
327     zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
328     end do
329     end do
330     end do
331     !
332     do jm=1,nd
333     do im=1,nd
334     do il=1,ncum
335     if(zm(il,im).ne.0.) then
336     ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
337     endif
338     end do
339     end do
340     end do
341     !
342     do jm=1,nd
343     do im=1,nd
344     do 999 il=1,ncum
345     qents(il,im,jm)=qent(il,im,jm)
346     ments(il,im,jm)=ment(il,im,jm)
347     999 continue
348     enddo
349     enddo
350    
351     return
352     end

  ViewVC Help
Powered by ViewVC 1.1.21