/[lmdze]/trunk/libf/phylmd/CV3_routines/cv3_mixing.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/CV3_routines/cv3_mixing.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (hide annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 11 months ago) by guez
File size: 11670 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

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     use cvparam3
7     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     integer ncum, nd, na, ntra, nloc
19     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     ! ori do 360 i=1,ncum*nlp
55     do 361 j=1,nl
56     do 360 i=1,ncum
57     nent(i,j)=0
58     ! in convect3, m is computed in cv3_closure
59     ! ori m(i,1)=0.0
60     360 continue
61     361 continue
62    
63     ! ori do 400 k=1,nlp
64     ! ori do 390 j=1,nlp
65     do 400 j=1,nl
66     do 390 k=1,nl
67     do 385 i=1,ncum
68     qent(i,k,j)=rr(i,j)
69     uent(i,k,j)=u(i,j)
70     vent(i,k,j)=v(i,j)
71     elij(i,k,j)=0.0
72     !ym ment(i,k,j)=0.0
73     !ym sij(i,k,j)=0.0
74     385 continue
75     390 continue
76     400 continue
77    
78     !ym
79     ment(1:ncum,1:nd,1:nd)=0.0
80     sij(1:ncum,1:nd,1:nd)=0.0
81    
82     ! do k=1,ntra
83     ! do j=1,nd ! instead nlp
84     ! do i=1,nd ! instead nlp
85     ! do il=1,ncum
86     ! traent(il,i,j,k)=tra(il,j,k)
87     ! enddo
88     ! enddo
89     ! enddo
90     ! enddo
91     zm(:,:)=0.
92    
93     !=====================================================================
94     ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
95     ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
96     ! --- FRACTION (sij)
97     !=====================================================================
98    
99     do 750 i=minorig+1, nl
100    
101     do 710 j=minorig,nl
102     do 700 il=1,ncum
103     if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
104     (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
105    
106     rti=rr(il,1)-ep(il,i)*clw(il,i)
107     bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
108     anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
109     denom=h(il,i)-hp(il,i)+(cpd-cpv)*(rr(il,i)-rti)*t(il,j)
110     dei=denom
111     if(abs(dei).lt.0.01)dei=0.01
112     sij(il,i,j)=anum/dei
113     sij(il,i,i)=1.0
114     altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
115     altem=altem/bf2
116     cwat=clw(il,j)*(1.-ep(il,j))
117     stemp=sij(il,i,j)
118     if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat) &
119     .and.j.gt.i)then
120     anum=anum-lv(il,j)*(rti-rs(il,j)-cwat*bf2)
121     denom=denom+lv(il,j)*(rr(il,i)-rti)
122     if(abs(denom).lt.0.01)denom=0.01
123     sij(il,i,j)=anum/denom
124     altem=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti-rs(il,j)
125     altem=altem-(bf2-1.)*cwat
126     end if
127     if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
128     qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti
129     uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))
130     vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))
131     !!!! do k=1,ntra
132     !!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
133     !!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
134     !!!! end do
135     elij(il,i,j)=altem
136     elij(il,i,j)=amax1(0.0,elij(il,i,j))
137     ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
138     nent(il,i)=nent(il,i)+1
139     end if
140     sij(il,i,j)=amax1(0.0,sij(il,i,j))
141     sij(il,i,j)=amin1(1.0,sij(il,i,j))
142     endif ! new
143     700 continue
144     710 continue
145    
146     ! do k=1,ntra
147     ! do j=minorig,nl
148     ! do il=1,ncum
149     ! if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
150     ! : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
151     ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
152     ! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
153     ! endif
154     ! enddo
155     ! enddo
156     ! enddo
157    
158     !
159     ! *** if no air can entrain at level i assume that updraft detrains ***
160     ! *** at that level and calculate detrained air flux and properties ***
161     !
162    
163     !@ do 170 i=icb(il),inb(il)
164    
165     do 740 il=1,ncum
166     if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
167     !@ if(nent(il,i).eq.0)then
168     ment(il,i,i)=m(il,i)
169     qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
170     uent(il,i,i)=u(il,nk(il))
171     vent(il,i,i)=v(il,nk(il))
172     elij(il,i,i)=clw(il,i)
173     !MAF sij(il,i,i)=1.0
174     sij(il,i,i)=0.0
175     end if
176     740 continue
177     750 continue
178    
179     ! do j=1,ntra
180     ! do i=minorig+1,nl
181     ! do il=1,ncum
182     ! if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
183     ! traent(il,i,i,j)=tra(il,nk(il),j)
184     ! endif
185     ! enddo
186     ! enddo
187     ! enddo
188    
189     do 100 j=minorig,nl
190     do 101 i=minorig,nl
191     do 102 il=1,ncum
192     if ((j.ge.(icb(il)-1)).and.(j.le.inb(il)) &
193     .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
194     sigij(il,i,j)=sij(il,i,j)
195     endif
196     102 continue
197     101 continue
198     100 continue
199     !@ enddo
200    
201     !@170 continue
202    
203     !=====================================================================
204     ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
205     ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
206     !=====================================================================
207    
208     !ym call zilch(asum,ncum*nd)
209     !ym call zilch(bsum,ncum*nd)
210     !ym call zilch(csum,ncum*nd)
211     call zilch(asum,nloc*nd)
212     call zilch(csum,nloc*nd)
213     call zilch(csum,nloc*nd)
214    
215     do il=1,ncum
216     lwork(il) = .FALSE.
217     enddo
218    
219     DO 789 i=minorig+1,nl
220    
221     num1=0
222     do il=1,ncum
223     if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
224     enddo
225     if (num1.le.0) goto 789
226    
227    
228     do 781 il=1,ncum
229     if ( i.ge.icb(il) .and. i.le.inb(il) ) then
230     lwork(il)=(nent(il,i).ne.0)
231     qp=rr(il,1)-ep(il,i)*clw(il,i)
232     anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i)) &
233     +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
234     denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp) &
235     +(cpd-cpv)*t(il,i)*(rr(il,i)-qp)
236     if(abs(denom).lt.0.01)denom=0.01
237     scrit(il)=anum/denom
238     alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
239     if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
240     smax(il)=0.0
241     asij(il)=0.0
242     endif
243     781 continue
244    
245     do 175 j=nl,minorig,-1
246    
247     num2=0
248     do il=1,ncum
249     if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
250     j.ge.(icb(il)-1) .and. j.le.inb(il) &
251     .and. lwork(il) ) num2=num2+1
252     enddo
253     if (num2.le.0) goto 175
254    
255     do 782 il=1,ncum
256     if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
257     j.ge.(icb(il)-1) .and. j.le.inb(il) &
258     .and. lwork(il) ) then
259    
260     if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
261     wgh=1.0
262     if(j.gt.i)then
263     sjmax=amax1(sij(il,i,j+1),smax(il))
264     sjmax=amin1(sjmax,scrit(il))
265     smax(il)=amax1(sij(il,i,j),smax(il))
266     sjmin=amax1(sij(il,i,j-1),smax(il))
267     sjmin=amin1(sjmin,scrit(il))
268     if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
269     smid=amin1(sij(il,i,j),scrit(il))
270     else
271     sjmax=amax1(sij(il,i,j+1),scrit(il))
272     smid=amax1(sij(il,i,j),scrit(il))
273     sjmin=0.0
274     if(j.gt.1)sjmin=sij(il,i,j-1)
275     sjmin=amax1(sjmin,scrit(il))
276     endif
277     delp=abs(sjmax-smid)
278     delm=abs(sjmin-smid)
279     asij(il)=asij(il)+wgh*(delp+delm)
280     ment(il,i,j)=ment(il,i,j)*(delp+delm)*wgh
281     endif
282     endif
283     782 continue
284    
285     175 continue
286    
287     do il=1,ncum
288     if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
289     asij(il)=amax1(1.0e-16,asij(il))
290     asij(il)=1.0/asij(il)
291     asum(il,i)=0.0
292     bsum(il,i)=0.0
293     csum(il,i)=0.0
294     endif
295     enddo
296    
297     do 180 j=minorig,nl
298     do il=1,ncum
299     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
300     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
301     ment(il,i,j)=ment(il,i,j)*asij(il)
302     endif
303     enddo
304     180 continue
305    
306     do 190 j=minorig,nl
307     do il=1,ncum
308     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
309     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
310     asum(il,i)=asum(il,i)+ment(il,i,j)
311     ment(il,i,j)=ment(il,i,j)*sig(il,j)
312     bsum(il,i)=bsum(il,i)+ment(il,i,j)
313     endif
314     enddo
315     190 continue
316    
317     do il=1,ncum
318     if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
319     bsum(il,i)=amax1(bsum(il,i),1.0e-16)
320     bsum(il,i)=1.0/bsum(il,i)
321     endif
322     enddo
323    
324     do 195 j=minorig,nl
325     do il=1,ncum
326     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
327     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
328     ment(il,i,j)=ment(il,i,j)*asum(il,i)*bsum(il,i)
329     endif
330     enddo
331     195 continue
332    
333     do 197 j=minorig,nl
334     do il=1,ncum
335     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
336     .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
337     csum(il,i)=csum(il,i)+ment(il,i,j)
338     endif
339     enddo
340     197 continue
341    
342     do il=1,ncum
343     if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il) &
344     .and. csum(il,i).lt.m(il,i) ) then
345     nent(il,i)=0
346     ment(il,i,i)=m(il,i)
347     qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
348     uent(il,i,i)=u(il,nk(il))
349     vent(il,i,i)=v(il,nk(il))
350     elij(il,i,i)=clw(il,i)
351     !MAF sij(il,i,i)=1.0
352     sij(il,i,i)=0.0
353     endif
354     enddo ! il
355    
356     ! do j=1,ntra
357     ! do il=1,ncum
358     ! if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
359     ! : .and. csum(il,i).lt.m(il,i) ) then
360     ! traent(il,i,i,j)=tra(il,nk(il),j)
361     ! endif
362     ! enddo
363     ! enddo
364     789 continue
365     !
366     ! MAF: renormalisation de MENT
367     do jm=1,nd
368     do im=1,nd
369     do il=1,ncum
370     zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
371     end do
372     end do
373     end do
374     !
375     do jm=1,nd
376     do im=1,nd
377     do il=1,ncum
378     if(zm(il,im).ne.0.) then
379     ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
380     endif
381     end do
382     end do
383     end do
384     !
385     do jm=1,nd
386     do im=1,nd
387     do 999 il=1,ncum
388     qents(il,im,jm)=qent(il,im,jm)
389     ments(il,im,jm)=ment(il,im,jm)
390     999 continue
391     enddo
392     enddo
393    
394     return
395     end

  ViewVC Help
Powered by ViewVC 1.1.21