/[lmdze]/trunk/libf/phylmd/CV_routines/cv_yield.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/CV_routines/cv_yield.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (show annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 7 months ago) by guez
File size: 13307 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

1
2 SUBROUTINE cv_yield(nloc,ncum,nd,nk,icb,inb,delt &
3 ,t,q,u,v,gz,p,ph,h,hp,lv,cpn &
4 ,ep,clw,frac,m,mp,qp,up,vp &
5 ,wt,water,evap &
6 ,ment,qent,uent,vent,nent,elij &
7 ,tv,tvp &
8 ,iflag,wd,qprime,tprime &
9 ,precip,cbmf,ft,fq,fu,fv,Ma,qcondc)
10 use cvthermo
11 use cvparam
12 implicit none
13
14
15 ! inputs
16 integer ncum, nd, nloc
17 integer nk(nloc), icb(nloc), inb(nloc)
18 integer nent(nloc,nd)
19 real, intent(in):: delt
20 real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd)
21 real gz(nloc,nd)
22 real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
23 real hp(nloc,nd), lv(nloc,nd)
24 real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc)
25 real m(nloc,nd), mp(nloc,nd), qp(nloc,nd)
26 real up(nloc,nd), vp(nloc,nd)
27 real wt(nloc,nd), water(nloc,nd), evap(nloc,nd)
28 real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd)
29 real uent(nloc,nd,nd), vent(nloc,nd,nd)
30 real tv(nloc,nd), tvp(nloc,nd)
31
32 ! outputs
33 integer iflag(nloc) ! also an input
34 real cbmf(nloc) ! also an input
35 real wd(nloc), tprime(nloc), qprime(nloc)
36 real precip(nloc)
37 real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
38 real Ma(nloc,nd)
39 real qcondc(nloc,nd)
40
41 ! local variables
42 integer i,j,ij,k,num1
43 real dpinv,cpinv,awat,fqold,ftold,fuold,fvold,delti
44 real work(nloc), am(nloc),amp1(nloc),ad(nloc)
45 real ents(nloc), uav(nloc),vav(nloc),lvcp(nloc,nd)
46 real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
47 real siga(nloc,nd), ax(nloc,nd), mac(nloc,nd) ! cld
48
49
50 ! -- initializations:
51
52 delti = 1.0/delt
53
54 do 160 i=1,ncum
55 precip(i)=0.0
56 wd(i)=0.0
57 tprime(i)=0.0
58 qprime(i)=0.0
59 do 170 k=1,nl+1
60 ft(i,k)=0.0
61 fu(i,k)=0.0
62 fv(i,k)=0.0
63 fq(i,k)=0.0
64 lvcp(i,k)=lv(i,k)/cpn(i,k)
65 qcondc(i,k)=0.0 ! cld
66 qcond(i,k)=0.0 ! cld
67 nqcond(i,k)=0.0 ! cld
68 170 continue
69 160 continue
70
71 !
72 ! *** Calculate surface precipitation in mm/day ***
73 !
74 do 1190 i=1,ncum
75 if(iflag(i).le.1)then
76 precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
77 endif
78 1190 continue
79 !
80 !
81 ! *** Calculate downdraft velocity scale and surface temperature and ***
82 ! *** water vapor fluctuations ***
83 !
84 do i=1,ncum
85 wd(i)=betad*abs(mp(i,icb(i)))*0.01*rrd*t(i,icb(i)) &
86 /(sigd*p(i,icb(i)))
87 qprime(i)=0.5*(qp(i,1)-q(i,1))
88 tprime(i)=lv0*qprime(i)/cpd
89 enddo
90 !
91 ! *** Calculate tendencies of lowest level potential temperature ***
92 ! *** and mixing ratio ***
93 !
94 do 1200 i=1,ncum
95 work(i)=0.01/(ph(i,1)-ph(i,2))
96 am(i)=0.0
97 1200 continue
98 do 1220 k=2,nl
99 do 1210 i=1,ncum
100 if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
101 am(i)=am(i)+m(i,k)
102 endif
103 1210 continue
104 1220 continue
105 do 1240 i=1,ncum
106 if((g*work(i)*am(i)).ge.delti)iflag(i)=1
107 ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1) &
108 +(gz(i,2)-gz(i,1))/cpn(i,1))
109 ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1)
110 ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2) &
111 -t(i,1))*work(i)/cpn(i,1)
112 fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))* &
113 work(i)+sigd*evap(i,1)
114 fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i)
115 fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1)) &
116 +am(i)*(u(i,2)-u(i,1)))
117 fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1)) &
118 +am(i)*(v(i,2)-v(i,1)))
119 1240 continue
120 do 1260 j=2,nl
121 do 1250 i=1,ncum
122 if(j.le.inb(i))then
123 fq(i,1)=fq(i,1) &
124 +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1))
125 fu(i,1)=fu(i,1) &
126 +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1))
127 fv(i,1)=fv(i,1) &
128 +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1))
129 endif
130 1250 continue
131 1260 continue
132 !
133 ! *** Calculate tendencies of potential temperature and mixing ratio ***
134 ! *** at levels above the lowest level ***
135 !
136 ! *** First find the net saturated updraft and downdraft mass fluxes ***
137 ! *** through each level ***
138 !
139 do 1500 i=2,nl+1
140 !
141 num1=0
142 do 1265 ij=1,ncum
143 if(i.le.inb(ij))num1=num1+1
144 1265 continue
145 if(num1.le.0)go to 1500
146 !
147 call zilch(amp1,ncum)
148 call zilch(ad,ncum)
149 !
150 do 1280 k=i+1,nl+1
151 do 1270 ij=1,ncum
152 if((i.ge.nk(ij)).and.(i.le.inb(ij)) &
153 .and.(k.le.(inb(ij)+1)))then
154 amp1(ij)=amp1(ij)+m(ij,k)
155 endif
156 1270 continue
157 1280 continue
158 !
159 do 1310 k=1,i
160 do 1300 j=i+1,nl+1
161 do 1290 ij=1,ncum
162 if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
163 amp1(ij)=amp1(ij)+ment(ij,k,j)
164 endif
165 1290 continue
166 1300 continue
167 1310 continue
168 do 1340 k=1,i-1
169 do 1330 j=i,nl+1
170 do 1320 ij=1,ncum
171 if((i.le.inb(ij)).and.(j.le.inb(ij)))then
172 ad(ij)=ad(ij)+ment(ij,j,k)
173 endif
174 1320 continue
175 1330 continue
176 1340 continue
177 !
178 do 1350 ij=1,ncum
179 if(i.le.inb(ij))then
180 dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
181 cpinv=1.0/cpn(ij,i)
182 !
183 ft(ij,i)=ft(ij,i) &
184 +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i) &
185 +(gz(ij,i+1)-gz(ij,i))*cpinv) &
186 -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) &
187 -sigd*lvcp(ij,i)*evap(ij,i)
188 ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+ &
189 t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
190 ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)* &
191 (t(ij,i+1)-t(ij,i))*dpinv*cpinv
192 fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))- &
193 ad(ij)*(q(ij,i)-q(ij,i-1)))
194 fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))- &
195 ad(ij)*(u(ij,i)-u(ij,i-1)))
196 fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))- &
197 ad(ij)*(v(ij,i)-v(ij,i-1)))
198 endif
199 1350 continue
200 do 1370 k=1,i-1
201 do 1360 ij=1,ncum
202 if(i.le.inb(ij))then
203 awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
204 awat=max(awat,0.0)
205 fq(ij,i)=fq(ij,i) &
206 +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
207 fu(ij,i)=fu(ij,i) &
208 +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
209 fv(ij,i)=fv(ij,i) &
210 +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
211 ! (saturated updrafts resulting from mixing) ! cld
212 qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld
213 nqcond(ij,i)=nqcond(ij,i)+1. ! cld
214 endif
215 1360 continue
216 1370 continue
217 do 1390 k=i,nl+1
218 do 1380 ij=1,ncum
219 if((i.le.inb(ij)).and.(k.le.inb(ij)))then
220 fq(ij,i)=fq(ij,i) &
221 +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
222 fu(ij,i)=fu(ij,i) &
223 +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
224 fv(ij,i)=fv(ij,i) &
225 +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
226 endif
227 1380 continue
228 1390 continue
229 do 1400 ij=1,ncum
230 if(i.le.inb(ij))then
231 fq(ij,i)=fq(ij,i) &
232 +sigd*evap(ij,i)+g*(mp(ij,i+1)* &
233 (qp(ij,i+1)-q(ij,i)) &
234 -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
235 fu(ij,i)=fu(ij,i) &
236 +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)* &
237 (up(ij,i)-u(ij,i-1)))*dpinv
238 fv(ij,i)=fv(ij,i) &
239 +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)* &
240 (vp(ij,i)-v(ij,i-1)))*dpinv
241 ! (saturated downdrafts resulting from mixing) ! cld
242 do k=i+1,inb(ij) ! cld
243 qcond(ij,i)=qcond(ij,i)+elij(ij,k,i) ! cld
244 nqcond(ij,i)=nqcond(ij,i)+1. ! cld
245 enddo ! cld
246 ! (particular case: no detraining level is found) ! cld
247 if (nent(ij,i).eq.0) then ! cld
248 qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld
249 nqcond(ij,i)=nqcond(ij,i)+1. ! cld
250 endif ! cld
251 if (nqcond(ij,i).ne.0.) then ! cld
252 qcond(ij,i)=qcond(ij,i)/nqcond(ij,i) ! cld
253 endif ! cld
254 endif
255 1400 continue
256 1500 continue
257 !
258 ! *** Adjust tendencies at top of convection layer to reflect ***
259 ! *** actual position of the level zero cape ***
260 !
261 do 503 ij=1,ncum
262 fqold=fq(ij,inb(ij))
263 fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
264 fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1) &
265 +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ &
266 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij)) &
267 /lv(ij,inb(ij)-1)
268 ftold=ft(ij,inb(ij))
269 ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
270 ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1) &
271 +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ &
272 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij)) &
273 /cpn(ij,inb(ij)-1)
274 fuold=fu(ij,inb(ij))
275 fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
276 fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1) &
277 +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ &
278 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
279 fvold=fv(ij,inb(ij))
280 fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
281 fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1) &
282 +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/ &
283 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
284 503 continue
285 !
286 ! *** Very slightly adjust tendencies to force exact ***
287 ! *** enthalpy, momentum and tracer conservation ***
288 !
289 do 682 ij=1,ncum
290 ents(ij)=0.0
291 uav(ij)=0.0
292 vav(ij)=0.0
293 do 681 i=1,inb(ij)
294 ents(ij)=ents(ij) &
295 +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1))
296 uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
297 vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
298 681 continue
299 682 continue
300 do 683 ij=1,ncum
301 ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
302 uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
303 vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
304 683 continue
305 do 642 ij=1,ncum
306 do 641 i=1,inb(ij)
307 ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
308 fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
309 fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
310 641 continue
311 642 continue
312 !
313 do 1810 k=1,nl+1
314 do 1800 i=1,ncum
315 if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
316 1800 continue
317 1810 continue
318 !
319 !
320 do 1900 i=1,ncum
321 if(iflag(i).gt.2)then
322 precip(i)=0.0
323 cbmf(i)=0.0
324 endif
325 1900 continue
326 do 1920 k=1,nl
327 do 1910 i=1,ncum
328 if(iflag(i).gt.2)then
329 ft(i,k)=0.0
330 fq(i,k)=0.0
331 fu(i,k)=0.0
332 fv(i,k)=0.0
333 qcondc(i,k)=0.0 ! cld
334 endif
335 1910 continue
336 1920 continue
337
338 do k=1,nl+1
339 do i=1,ncum
340 Ma(i,k) = 0.
341 enddo
342 enddo
343 do k=nl,1,-1
344 do i=1,ncum
345 Ma(i,k) = Ma(i,k+1)+m(i,k)
346 enddo
347 enddo
348
349 !
350 ! *** diagnose the in-cloud mixing ratio ***
351 ! *** of condensed water ***
352 !
353 DO ij=1,ncum
354 do i=1,nd
355 mac(ij,i)=0.0
356 wa(ij,i)=0.0
357 siga(ij,i)=0.0
358 enddo
359 do i=nk(ij),inb(ij)
360 do k=i+1,inb(ij)+1
361 mac(ij,i)=mac(ij,i)+m(ij,k)
362 enddo
363 enddo
364 do i=icb(ij),inb(ij)-1
365 ax(ij,i)=0.
366 do j=icb(ij),i
367 ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j)) &
368 *(ph(ij,j)-ph(ij,j+1))/p(ij,j)
369 enddo
370 if (ax(ij,i).gt.0.0) then
371 wa(ij,i)=sqrt(2.*ax(ij,i))
372 endif
373 enddo
374 do i=1,nl
375 if (wa(ij,i).gt.0.0) &
376 siga(ij,i)=mac(ij,i)/wa(ij,i) &
377 *rrd*tvp(ij,i)/p(ij,i)/100./delta
378 siga(ij,i) = min(siga(ij,i),1.0)
379 qcondc(ij,i)=siga(ij,i)*clw(ij,i)*(1.-ep(ij,i)) &
380 + (1.-siga(ij,i))*qcond(ij,i)
381 enddo
382 ENDDO
383
384 return
385 end

  ViewVC Help
Powered by ViewVC 1.1.21