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

  ViewVC Help
Powered by ViewVC 1.1.21