/[lmdze]/trunk/phylmd/CV_routines/cv_mixing.f
ViewVC logotype

Annotation of /trunk/phylmd/CV_routines/cv_mixing.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 103 - (hide annotations)
Fri Aug 29 13:00:05 2014 UTC (9 years, 9 months ago) by guez
File size: 9408 byte(s) 
Renamed module cvparam to cv_param. Deleted procedure
cv_param. Changed variables of module cv_param into parameters.

In procedures cv_driver, cv_uncompress and cv3_uncompress, removed
some arguments giving dimensions and used module variables klon and
klev instead.

In procedures gradiv2, laplacien_gam and laplacien, changed
declarations of local variables because klevel is not always klev.

Removed code for nudging surface pressure.

Removed arguments pim and pjm of tau2alpha. Added assignment of false
to variable first.

Replaced real argument del of procedures foeew and FOEDE by logical
argument.
1 guez 52
2     SUBROUTINE cv_mixing(nloc,ncum,nd,icb,nk,inb,inb1 &
3     ,ph,t,q,qs,u,v,h,lv,qnk &
4     ,hp,tv,tvp,ep,clw,cbmf &
5     ,m,ment,qent,uent,vent,nent,sij,elij)
6     use cvthermo
7 guez 103 use cv_param
8 guez 52 implicit none
9    
10    
11     ! inputs:
12 guez 97 integer, intent(in):: ncum, nd, nloc
13 guez 52 integer icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
14     real cbmf(nloc), qnk(nloc)
15     real ph(nloc,nd+1)
16     real t(nloc,nd), q(nloc,nd), qs(nloc,nd), lv(nloc,nd)
17     real u(nloc,nd), v(nloc,nd), h(nloc,nd), hp(nloc,nd)
18     real tv(nloc,nd), tvp(nloc,nd), ep(nloc,nd), clw(nloc,nd)
19    
20     ! outputs:
21     integer nent(nloc,nd)
22     real m(nloc,nd), ment(nloc,nd,nd), qent(nloc,nd,nd)
23     real uent(nloc,nd,nd), vent(nloc,nd,nd)
24     real sij(nloc,nd,nd), elij(nloc,nd,nd)
25    
26     ! local variables:
27     integer i, j, k, ij
28     integer num1, num2
29     real dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
30     real alt, qp1, smid, sjmin, sjmax, delp, delm
31     real work(nloc), asij(nloc), smin(nloc), scrit(nloc)
32     real bsum(nloc,nd)
33     logical lwork(nloc)
34    
35     !=====================================================================
36     ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
37     !=====================================================================
38     !
39     do 360 i=1,ncum*nlp
40     nent(i,1)=0
41     m(i,1)=0.0
42     360 continue
43     !
44     do 400 k=1,nlp
45     do 390 j=1,nlp
46     do 385 i=1,ncum
47     qent(i,k,j)=q(i,j)
48     uent(i,k,j)=u(i,j)
49     vent(i,k,j)=v(i,j)
50     elij(i,k,j)=0.0
51     ment(i,k,j)=0.0
52     sij(i,k,j)=0.0
53     385 continue
54     390 continue
55     400 continue
56     !
57     !-------------------------------------------------------------------
58     ! --- Calculate rates of mixing, m(i)
59     !-------------------------------------------------------------------
60     !
61     call zilch(work,ncum)
62     !
63     do 670 j=minorig+1,nl
64     do 660 i=1,ncum
65     if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
66     k=min(j,inb1(i))
67     dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) &
68     +entp*0.04*(ph(i,k)-ph(i,k+1))
69     work(i)=work(i)+dbo
70     m(i,j)=cbmf(i)*dbo
71     endif
72     660 continue
73     670 continue
74     do 690 k=minorig+1,nl
75     do 680 i=1,ncum
76     if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
77     m(i,k)=m(i,k)/work(i)
78     endif
79     680 continue
80     690 continue
81     !
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     !
90     do 750 i=minorig+1,nl
91     do 710 j=minorig+1,nl
92     do 700 ij=1,ncum
93     if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij)) &
94     .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
95     qti=qnk(ij)-ep(ij,i)*clw(ij,i)
96     bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j) &
97     /(rrv*t(ij,j)*t(ij,j)*cpd)
98     anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j))
99     denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j)
100     dei=denom
101     if(abs(dei).lt.0.01)dei=0.01
102     sij(ij,i,j)=anum/dei
103     sij(ij,i,i)=1.0
104     altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
105     altem=altem/bf2
106     cwat=clw(ij,j)*(1.-ep(ij,j))
107     stemp=sij(ij,i,j)
108     if((stemp.lt.0.0.or.stemp.gt.1.0.or. &
109     altem.gt.cwat).and.j.gt.i)then
110     anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2)
111     denom=denom+lv(ij,j)*(q(ij,i)-qti)
112     if(abs(denom).lt.0.01)denom=0.01
113     sij(ij,i,j)=anum/denom
114     altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
115     altem=altem-(bf2-1.)*cwat
116     endif
117     if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
118     qent(ij,i,j)=sij(ij,i,j)*q(ij,i) &
119     +(1.-sij(ij,i,j))*qti
120     uent(ij,i,j)=sij(ij,i,j)*u(ij,i) &
121     +(1.-sij(ij,i,j))*u(ij,nk(ij))
122     vent(ij,i,j)=sij(ij,i,j)*v(ij,i) &
123     +(1.-sij(ij,i,j))*v(ij,nk(ij))
124     elij(ij,i,j)=altem
125     elij(ij,i,j)=max(0.0,elij(ij,i,j))
126     ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
127     nent(ij,i)=nent(ij,i)+1
128     endif
129     sij(ij,i,j)=max(0.0,sij(ij,i,j))
130     sij(ij,i,j)=min(1.0,sij(ij,i,j))
131     endif
132     700 continue
133     710 continue
134     !
135     ! *** If no air can entrain at level i assume that updraft detrains ***
136     ! *** at that level and calculate detrained air flux and properties ***
137     !
138     do 740 ij=1,ncum
139     if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij)) &
140     .and.(nent(ij,i).eq.0))then
141     ment(ij,i,i)=m(ij,i)
142     qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
143     uent(ij,i,i)=u(ij,nk(ij))
144     vent(ij,i,i)=v(ij,nk(ij))
145     elij(ij,i,i)=clw(ij,i)
146     sij(ij,i,i)=1.0
147     endif
148     740 continue
149     750 continue
150     !
151     do 770 i=1,ncum
152     sij(i,inb(i),inb(i))=1.0
153     770 continue
154     !
155     !=====================================================================
156     ! --- NORMALIZE ENTRAINED AIR MASS FLUXES
157     ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
158     !=====================================================================
159     !
160     call zilch(bsum,ncum*nlp)
161     do 780 ij=1,ncum
162     lwork(ij)=.false.
163     780 continue
164     do 789 i=minorig+1,nl
165     !
166     num1=0
167     do 779 ij=1,ncum
168     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
169     779 continue
170     if(num1.le.0)go to 789
171     !
172     do 781 ij=1,ncum
173     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
174     lwork(ij)=(nent(ij,i).ne.0)
175     qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
176     anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
177     denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
178     if(abs(denom).lt.0.01)denom=0.01
179     scrit(ij)=anum/denom
180     alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
181     if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
182     asij(ij)=0.0
183     smin(ij)=1.0
184     endif
185     781 continue
186     do 783 j=minorig,nl
187     !
188     num2=0
189     do 778 ij=1,ncum
190     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) &
191     .and.(j.ge.icb(ij)).and.(j.le.inb(ij)) &
192     .and.lwork(ij))num2=num2+1
193     778 continue
194     if(num2.le.0)go to 783
195     !
196     do 782 ij=1,ncum
197     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) &
198     .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
199     if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
200     if(j.gt.i)then
201     smid=min(sij(ij,i,j),scrit(ij))
202     sjmax=smid
203     sjmin=smid
204     if(smid.lt.smin(ij) &
205     .and.sij(ij,i,j+1).lt.smid)then
206     smin(ij)=smid
207     sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
208     sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
209     sjmin=min(sjmin,scrit(ij))
210     endif
211     else
212     sjmax=max(sij(ij,i,j+1),scrit(ij))
213     smid=max(sij(ij,i,j),scrit(ij))
214     sjmin=0.0
215     if(j.gt.1)sjmin=sij(ij,i,j-1)
216     sjmin=max(sjmin,scrit(ij))
217     endif
218     delp=abs(sjmax-smid)
219     delm=abs(sjmin-smid)
220     asij(ij)=asij(ij)+(delp+delm) &
221     *(ph(ij,j)-ph(ij,j+1))
222     ment(ij,i,j)=ment(ij,i,j)*(delp+delm) &
223     *(ph(ij,j)-ph(ij,j+1))
224     endif
225     endif
226     782 continue
227     783 continue
228     do 784 ij=1,ncum
229     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
230     asij(ij)=max(1.0e-21,asij(ij))
231     asij(ij)=1.0/asij(ij)
232     bsum(ij,i)=0.0
233     endif
234     784 continue
235     do 786 j=minorig,nl+1
236     do 785 ij=1,ncum
237     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) &
238     .and.(j.ge.icb(ij)).and.(j.le.inb(ij)) &
239     .and.lwork(ij))then
240     ment(ij,i,j)=ment(ij,i,j)*asij(ij)
241     bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
242     endif
243     785 continue
244     786 continue
245     do 787 ij=1,ncum
246     if((i.ge.icb(ij)+1).and.(i.le.inb(ij)) &
247     .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
248     nent(ij,i)=0
249     ment(ij,i,i)=m(ij,i)
250     qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
251     uent(ij,i,i)=u(ij,nk(ij))
252     vent(ij,i,i)=v(ij,nk(ij))
253     elij(ij,i,i)=clw(ij,i)
254     sij(ij,i,i)=1.0
255     endif
256     787 continue
257     789 continue
258     !
259     return
260     end
  ViewVC Help
Powered by ViewVC 1.1.21