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

Contents of /trunk/phylmd/CV_routines/cv_yield.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21