Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

convmix_kfeta.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 22.0 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23	subroutine convmix_kfeta(itime)
24	! i
25	!**************************************************************
26	! handles all the calculations related to convective mixing
27	! Petra Seibert, Bernd C. Krueger, Feb 2001
28	! nested grids included, Bernd C. Krueger, May 2001
29	!
30	! Changes by Caroline Forster, April 2004 - February 2005:
31	! convmix called every lsynctime seconds
32	! CHANGES by A. Stohl:
33	! various run-time optimizations - February 2005
34
35	! CHANGED on 10.10.2007 save convective mass fluxes, update them every dt_conv
36
37	! CHANGED by Weiguo WANG 13 Aug, 2007, use KFeta CU convection scheme
38	!
39	! Changes by J. Brioude: particles sorting is much more efficient.
40	!
41	! input for kftea cumulus scheme
42	! u1d - 1-d u wind velocity profile
43	! v1d - 1-d v wind velocity profile
44	! t1d - 1-D temperture (K)
45	! qv1d- 1-D water vapor mixin gratio (kg/kg)
46	! p1d - 1-D pressure profile (pa)
47	! rho1d-1-D density profile (kg/m3)
48	! w0avg1d - 1-D vertical velocity (m/s)
49	! all above are defined at T-point or P-poit
50	! dz1d - dz between full levels
51	! delx - grid size of column (m)
52	! dt - integraiton time step (s)
53	! cudt - cumulus activation time interval (min)
54	! kts - starting point in z for convection calculation
55	! kte - ending point
56
57	! output
58	! umf - updraft mass flux
59	! uer - updraft entrainment flux
60	! udr - updraft detrainment flux
61	! dmf - downdraft mass flux
62	! der - downdraft entrainemnt flux
63	! ddr - downdraft detrainment flux
64	! cu_top1 -- top of cumulus cloud (index? for height)
65	! cu_bot1 -- bottom of cumulus cloud (index)
66
67	!**************************************************************
68
69	use flux_mod
70	use par_mod
71	use com_mod
72	use conv_mod
73
74	implicit none
75
76	integer :: igr,igrold, ipart, itime, ix, j, inest
77	integer :: ipconv
78	integer :: jy, kpart, ktop, ngrid,kz,kzp,a
79	integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
80	! itime [s] current time
81	! igrid(maxpart) horizontal grid position of each particle
82	! igridn(maxpart,maxnests) dto. for nested grids
83	! ipoint(maxpart) pointer to access particles according to grid position
84
85	logical :: lconv,warm_rain
86	real :: x, y, xtn,ytn, ztold, delt
87	real :: dt1,dt2,dtt,dummy
88	real :: duma, dumz(nuvzmax+1)
89	integer :: mind1,mind2
90	! dt1,dt2,dtt,mind1,mind2 variables used for time interpolation
91	integer :: itage,nage,duminc
92	real,parameter :: eps=nxmax/3.e5
93
94	integer :: i,k,numberp(maxpart)
95	!-- for KFeta
96	real, dimension(nuvzmax):: u1d,v1d,t1d,dz1d,qv1d,p1d, &
97	rho1d,w0avg1d,umf,uer,udr, &
98	dmf,der,ddr,zh
99	real :: cudt,delx,dt,cu_bot1,cu_top1,zp,fmix
100	real, dimension(nuvzmax+1):: umfzf,dmfzf,zf
101	integer :: kts,kte,if_update
102	!
103
104	! monitoring variables
105	real :: sumconv,sumall,sumpart
106	integer :: sumpartgrid(1000000)
107
108	! print *, "IN convmix_kfeta"
109	! write(*,'(//a,a//)')
110	! & '* Stopping in subr. convmix *',
111	! & ' This is not implemented for FLEXPART_WRF'
112	! stop
113
114
115	! Calculate auxiliary variables for time interpolation
116	!*****************************************************
117
118	delt=real(abs(lsynctime))
119
120	! dt_conv is given from input namelist
121	! dt_conv=3600.0
122	if_update=0
123	if ( mod(real(itime),dt_conv) .eq. 0 ) if_update=1
124	! print*,'conv itime',itime,dt_conv,if_update
125
126	! delt=dt_conv
127
128
129	dt1=real(itime-memtime(1))
130	dt2=real(memtime(2)-itime)
131	dtt=1./(dt1+dt2)
132	mind1=memind(1)
133	mind2=memind(2)
134
135	lconv = .false.
136
137	! for KFeta
138	! warm_rain=.true. ! depends on mp_physics in WRF, may add an option in the future
139	cudt = 10.0 ! cumulus para is called every 10 min in a time step, if dt < cudt, call once
140
141	kts=1
142	kte=nuvz-1
143
144	! if no particles are present return after initialization
145	!********************************************************
146
147	if (numpart.le.0) return
148
149	! Assign igrid and igridn, which are pseudo grid numbers indicating particles
150	! that are outside the part of the grid under consideration
151	! (e.g. particles near the poles or particles in other nests).
152	! Do this for all nests but use only the innermost nest; for all others
153	! igrid shall be -1
154	! Also, initialize index vector ipoint
155	!************************************************************************
156	! print*,'step1'
157	do ipart=1,numpart
158	igrid(ipart)=-1
159	do j=numbnests,1,-1
160	igridn(ipart,j)=-1
161	enddo
162	ipoint(ipart)=ipart
163	! do not consider particles that are (yet) not part of simulation
164	if (itra1(ipart).ne.itime) goto 20
165	x = xtra1(ipart)
166	y = ytra1(ipart)
167
168	! Determine which nesting level to be used
169	!**********************************************************
170
171	ngrid=0
172	do j=numbnests,1,-1
173	if ( x.gt.xln(j) .and. x.lt.xrn(j) .and. &
174	y.gt.yln(j) .and. y.lt.yrn(j) ) then
175	ngrid=j
176	goto 23
177	endif
178	end do
179	23 continue
180
181	! Determine nested grid coordinates
182	!**********************************
183
184	if (ngrid.gt.0) then
185	! nested grids
186	xtn=(x-xln(ngrid))*xresoln(ngrid)
187	ytn=(y-yln(ngrid))*yresoln(ngrid)
188	ix=nint(xtn)
189	jy=nint(ytn)
190	igridn(ipart,ngrid) = 1 + jy*nxn(ngrid) + ix
191	else if(ngrid.eq.0) then
192	! mother grid
193	ix=nint(x)
194	jy=nint(y)
195	igrid(ipart) = 1 + jy*nx + ix
196	endif
197
198	20 continue
199	end do
200
201	sumpart = 0.
202	sumconv = 0.
203
204	! print*,'step2'
205
206	!**************************************************************************************
207	! 1. Now, do everything for the mother domain and, later, for all of the nested domains
208	! While all particles have to be considered for redistribution, the Emanuel convection
209	! scheme only needs to be called once for every grid column where particles are present.
210	! Therefore, particles are sorted according to their grid position. Whenever a new grid
211	! cell is encountered by looping through the sorted particles, the convection scheme is called.
212	!**************************************************************************************
213	delx = dx
214
215	! sort particles according to horizontal position and calculate index vector IPOINT
216
217	call sort2(numpart,igrid,ipoint)
218	! print*,'step2 after sort',minval(igrid),maxval(igrid),numpart
219
220	! count particle # in each column
221	! do 40 i=abs(minval(igrid)),maxval(igrid)
222	! sumpart=0.
223	! do 41 k=1,numpart
224	! if(igrid(k) .eq. i) then
225	! sumpart=sumpart+1
226	! endif
227	!41 continue
228	! do 42 k=1,numparc Changes by J. Brioude: the sort of particles is much more efficient.t
229	!42 if(igrid(k) .eq. i) numberp(k)=int(sumpart)
230	!40 continue
231	! JB
232	if (maxval(igrid).gt.1000000) then
233	print*,'too much x and y grid. modify convmix_kfeta.f'
234	stop
235	endif
236	do k=1,1000000
237	sumpartgrid(k)=0
238	enddo
239	do k=1,numpart
240	if (igrid(k).gt.0) sumpartgrid(igrid(k))=sumpartgrid(igrid(k))+1
241	enddo
242	do k=1,numpart
243	if (igrid(k).gt.0) then
244	numberp(k)=sumpartgrid(igrid(k))
245	else
246	numberp(k)=0
247	endif
248	enddo
249
250	! print*,'step3',numpart
251
252	! Now visit all grid columns where particles are present
253	! by going through the sorted particles
254
255	igrold = -1
256	a=0
257	do kpart=1,numpart
258	igr = igrid(kpart)
259	if (igr .eq. -1 .or. numberp(kpart).le.20 &
260	! if (igr .eq. -1
261	) goto 50
262	ipart = ipoint(kpart)
263
264	! sumall = sumall + 1
265	!c For one column, we only need to compute 1D met once
266
267	if (igr .ne. igrold) then
268	sumconv=sumconv+1
269	! we are in a new grid column
270	jy = (igr-1)/nx
271	ix = igr - jy*nx - 1
272	a=a+1
273	! print*,'a',a
274	! Interpolate all meteorological data needed for the convection scheme
275
276	do kz=1,nuvz-1 ! nconvlev+1
277	! FLEXPART_WRF - used add_sfc_level for the shifting
278	! for W, it is not shifted, make sure w is 'true' vertical velocity!
279
280	kzp = kz + add_sfc_level
281	u1d(kz)=(u_wrf(ix,jy,kzp,mind1)*dt2+ &
282	u_wrf(ix,jy,kzp,mind2)dt1)dtt
283	v1d(kz)=(v_wrf(ix,jy,kzp,mind1)*dt2+ &
284	v_wrf(ix,jy,kzp,mind2)dt1)dtt
285	t1d(kz)=(tth(ix,jy,kzp,mind1)*dt2+ &
286	tth(ix,jy,kzp,mind2)dt1)dtt
287	qv1d(kz)=(qvh(ix,jy,kzp,mind1)*dt2+ &
288	qvh(ix,jy,kzp,mind2)dt1)dtt
289	p1d(kz)=(pph(ix,jy,kzp,mind1)*dt2+ &
290	pph(ix,jy,kzp,mind2)dt1)dtt
291	dz1d(kz)=(zzh(ix,jy,kzp+1,mind1)*dt2+ &
292	zzh(ix,jy,kzp+1,mind2)dt1)dtt- &
293	(zzh(ix,jy,kzp,mind1)*dt2+ &
294	zzh(ix,jy,kzp,mind2)dt1)dtt
295	w0avg1d(kz)=(w_wrf(ix,jy,kz,mind1)*dt2+ &
296	w_wrf(ix,jy,kz,mind2)dt1)dtt+ &
297	(w_wrf(ix,jy,kz+1,mind1)*dt2+ &
298	w_wrf(ix,jy,kz+1,mind2)dt1)dtt
299	w0avg1d(kz)=0.5*w0avg1d(kz)
300	rho1d(kz)=p1d(kz)/ &
301	(t1d(kz)(1.0+0.608qv1d(kz))) &
302	/287.0
303
304	! write(*,'(1x,I10,10F10.2)')kz,u1d(kz),v1d(kz),w0avg1d(kz),
305	! & t1d(kz),qv1d(kz),p1d(kz)/100,dz1d(kz)
306
307
308	enddo
309
310	! -- old convective mass fluxes
311	do k=kts,kte
312	umf(k)=umf3(ix,jy,k)
313	uer(k)=uer3(ix,jy,k)
314	udr(k)=udr3(ix,jy,k)
315	dmf(k)=dmf3(ix,jy,k)
316	der(k)=der3(ix,jy,k)
317	ddr(k)=ddr3(ix,jy,k)
318	enddo
319	cu_top1=cu_top(ix,jy)
320	cu_bot1=cu_bot(ix,jy)
321
322
323	! write(,)'1-D wind'
324
325	! Calculate convection flux, updrought flux, entrainment, detrainment flux
326	! downdraft flux, entrainment ,detrainment flux
327	warm_rain=.false.
328	if (mp_physics .eq. 1) warm_rain = .true.
329	! -- Update fluxes
330	if (if_update .eq. 1 ) then !! if_update
331	! write(,)'update convective fluxes, itime=',itime/3600.
332	! print*,u1d(4:8),v1d(4:8),t1d(4:8),dz1d(4:8),qv1d(4:8)
333	! print*,p1d(4:8),rho1d(4:8),w0avg1d(4:8)
334	! print*,cudt,delx,dt_conv,warm_rain
335	! print*,'attend'
336	! pause
337	CALL KF_ETA(nuvzmax,u1d,v1d,t1d,dz1d,qv1d,p1d, & ! IN
338	rho1d,w0avg1d,cudt,delx,dt_conv,warm_rain,kts,kte, & ! IN
339	umf,uer,udr,dmf,der,ddr,cu_bot1,cu_top1) ! OUT
340	dummy=0.
341	do k=kts,kte
342	umf3(ix,jy,k)=umf(k)
343	uer3(ix,jy,k)=uer(k)
344	udr3(ix,jy,k)=udr(k)
345	dmf3(ix,jy,k)=dmf(k)
346	der3(ix,jy,k)=der(k)
347	ddr3(ix,jy,k)=ddr(k)
348	dummy=dummy+umf(k)+uer(k)+udr(k)
349	enddo
350	if (dummy.gt.0.) then
351	! print*,'dummy',dummy
352	duminc=1
353	else
354	duminc=0
355	endif
356	cu_top(ix,jy)=cu_top1
357	cu_bot(ix,jy)=cu_bot1
358	! if (cu_top1.gt.1.) print*,'cu h',cu_bot1,cu_top1
359	endif !! if_update
360	! if (a.gt.2000) then
361	! print*,'after kf_eta',a
362	! a=0
363	! endif
364
365	! write(,)'ix,jy=',ix,jy,itime
366	! write(,)'previous column part#=',sumpart
367	! write(,)'FLUX,k,umf,uer,udr,dmf,der,ddr'
368	! write(,)'cu_bot1,cu_top1=',cu_bot1,cu_top1
369	! if (cu_top1 .lt. cu_bot1) write(,)'umf=', umf(1),umf(10)
370
371	! do kz=kts,kte
372	! write(*,'(1x,I10,10E10.2)')kz,umf(kz),uer(kz),udr(kz),dmf(kz)
373	! & ,der(kz),ddr(kz),p1d(kz)/100
374
375	! enddo
376
377	sumpart=0
378	IF (cu_top1 .gt. cu_bot1+1 ) then ! lconv
379	! write(250+lconvection,*)'-1',itime,ix,jy
380	lconv = .true.
381
382	! Prepare data for redistributing particle
383
384	CALL pre_redist_kf(nuvzmax,nuvz,umf,dmf,dz1d,p1d,delx,delt, & ! IN
385	cu_bot1,cu_top1, & ! IN
386	zf,zh, & ! OUT
387	umfzf,dmfzf,fmix) ! OUT
388
389	else
390	lconv= .false.
391
392	ENDIF ! lconv
393
394
395	igrold = igr
396	ktop = 0
397	endif
398
399	sumpart=sumpart+1
400	! treat particle only if column has convection
401	if (lconv .eqv. .true.) then
402	! assign new vertical position to particle
403	! ztold=ztra1(ipart)
404	zp=ztra1(ipart)
405	!C write(,)'befrore convection zp= ',zp
406
407	! write(,)'part No =',sumpart
408	CALL redist_kf(lconvection,ldirect,delt,delx, & ! IN
409	dz1d,nzmax, nz,umf,uer,udr,dmf, & ! IN
410	der,ddr,rho1d, & ! IN
411	zf,zh, & ! IN
412	umfzf,dmfzf,fmix, & ! IN
413	zp) ! IN/OUT
414
415	if (zp .lt. 0.0) zp=-1.0*zp
416	if (zp .gt. height(nz)-0.5) &
417	zp = height(nz)-0.5
418	! if (abs(zp-ztra1(ipart)) .ge. 1e-5)
419	! &write(250+lconvection,*)ztra1(ipart),zp,zp-ztra1(ipart)
420	! if (duminc.eq.1)
421	! + print*,'true conv',dummy,zp-ztra1(ipart),cu_top1-cu_bot1
422
423	ztra1(ipart) = zp
424	!C write(,)'after convection, zp=',zp
425
426	!C OLD call redist(ipart,ktop,ipconv)
427	! if (ipconv.le.0) sumconv = sumconv+1
428
429	! Calculate the gross fluxes across layer interfaces
430	!***************************************************
431
432	if (iflux.eq.1) then
433	itage=abs(itra1(ipart)-itramem(ipart))
434	do nage=1,nageclass
435	if (itage.lt.lage(nage)) goto 37
436	enddo
437	37 continue
438	! print*,'step 4'
439
440	if (nage.le.nageclass) &
441	call calcfluxes(nage,ipart,real(xtra1(ipart)), &
442	real(ytra1(ipart)),ztold)
443	endif
444
445	endif !(lconv .eqv. .true)
446	enddo
447	50 continue
448
449	! write(,)'total convective columns=',sumconv,
450	! & 'time=', 1.0*itime/3600
451
452	!***********************************************************************************
453	! 2. Nested domains
454	!***********************************************************************************
455
456	! sort particles according to horizontal position and calculate index vector IPOINT
457
458	do inest=1,numbnests
459	delx = dxn(inest)
460	if (delx .le. 10000.0) goto 70 ! for small grid size, no need to do convection
461	do ipart=1,numpart
462	ipoint(ipart)=ipart
463	igrid(ipart) = igridn(ipart,inest)
464	enddo
465	call sort2(numpart,igrid,ipoint)
466
467	! print*,'step in nest'
468	! Now visit all grid columns where particles are present
469	! by going through the sorted particles
470
471	igrold = -1
472	do kpart=1,numpart
473	igr = igrid(kpart)
474	if (igr .eq. -1) goto 60
475	ipart = ipoint(kpart)
476	! sumall = sumall + 1
477	if (igr .ne. igrold) then
478	! we are in a new grid column
479	jy = (igr-1)/nxn(inest)
480	ix = igr - jy*nxn(inest) - 1
481
482	! Interpolate all meteorological data needed for the convection scheme
483
484	! Interpolate all meteorological data needed for the convection scheme
485
486	do kz=1,nuvz ! nconvlev+1
487	! FLEXPART_WRF - used add_sfc_level for the shifting
488	! for W, it is not shifted, make sure w is 'true' vertical velocity!
489
490	kzp = kz + add_sfc_level
491	u1d(kz)=(un_wrf(ix,jy,kzp,mind1,inest)*dt2+ &
492	un_wrf(ix,jy,kzp,mind2,inest)dt1)dtt
493	v1d(kz)=(vn_wrf(ix,jy,kzp,mind1,inest)*dt2+ &
494	vn_wrf(ix,jy,kzp,mind2,inest)dt1)dtt
495	t1d(kz)=(tthn(ix,jy,kzp,mind1,inest)*dt2+ &
496	tthn(ix,jy,kzp,mind2,inest)dt1)dtt
497	qv1d(kz)=(qvhn(ix,jy,kzp,mind1,inest)*dt2+ &
498	qvhn(ix,jy,kzp,mind2,inest)dt1)dtt
499	p1d(kz)=(pphn(ix,jy,kzp,mind1,inest)*dt2+ &
500	pphn(ix,jy,kzp,mind2,inest)dt1)dtt
501	dz1d(kz)=(zzhn(ix,jy,kzp+1,mind1,inest)*dt2+ &
502	zzhn(ix,jy,kzp+1,mind2,inest)dt1)dtt- &
503	(zzhn(ix,jy,kzp,mind1,inest)*dt2+ &
504	zzhn(ix,jy,kzp,mind2,inest)dt1)dtt
505	w0avg1d(kz)=(wn_wrf(ix,jy,kz,mind1,inest)*dt2+ &
506	wn_wrf(ix,jy,kz,mind2,inest)dt1)dtt+ &
507	(wn_wrf(ix,jy,kz+1,mind1,inest)*dt2+ &
508	wn_wrf(ix,jy,kz+1,mind2,inest)dt1)dtt
509	w0avg1d(kz)=0.5*w0avg1d(kz)
510	rho1d(kz)=p1d(kz)/ &
511	(t1d(kz)(1.0+0.608qv1d(kz))) &
512	/287.0
513
514	enddo
515
516	!C Old convective mass fluxes
517	do k=kts,kte
518	umf(k)=umf3n(ix,jy,k,inest)
519	uer(k)=uer3n(ix,jy,k,inest)
520	udr(k)=udr3n(ix,jy,k,inest)
521	dmf(k)=dmf3n(ix,jy,k,inest)
522	der(k)=der3n(ix,jy,k,inest)
523	ddr(k)=ddr3n(ix,jy,k,inest)
524	enddo
525	cu_top1=cu_topn(ix,jy,inest)
526	cu_bot1=cu_botn(ix,jy,inest)
527
528
529
530	!alculate convection flux, updrought flux, entrainment, detrainment flux
531	! downdraft flux, entrainment ,detrainment flux
532	warm_rain = .false.
533	if (mp_physicsn(inest) .eq. 1 ) warm_rain = .true.
534
535	if (if_update .eq. 1 ) then !!! update
536	CALL KF_ETA(nuvzmax,u1d,v1d,t1d,dz1d,qv1d,p1d, & ! IN
537	rho1d,w0avg1d,cudt,delx,dt_conv,warm_rain,kts,kte, & ! IN
538	umf,uer,udr,dmf,der,ddr,cu_bot1,cu_top1) ! OUT
539	do k=kts,kte
540	umf3n(ix,jy,k,inest)=umf(k)
541	uer3n(ix,jy,k,inest)=uer(k)
542	udr3n(ix,jy,k,inest)=udr(k)
543	dmf3n(ix,jy,k,inest)=dmf(k)
544	der3n(ix,jy,k,inest)=der(k)
545	ddr3n(ix,jy,k,inest)=ddr(k)
546	enddo
547	cu_topn(ix,jy,inest)=cu_top1
548	cu_botn(ix,jy,inest)=cu_bot1
549	endif !!! update
550
551
552
553	IF (cu_top1 .gt. cu_bot1) then ! lconv
554
555	lconv = .true.
556
557	! Prepare data for redistributing particle
558
559	CALL pre_redist_kf(nuvzmax,nuvz,umf,dmf,dz1d,p1d,delx,delt, & ! IN
560	cu_bot1,cu_top1, & ! IN
561	zf,zh, & ! OUT
562	umfzf,dmfzf,fmix) ! OUT
563
564	else
565	lconv = .false.
566	ENDIF ! lconv
567
568
569	igrold = igr
570	ktop = 0
571	endif
572
573	! treat particle only if column has convection
574	if (lconv .eqv. .true.) then
575	! assign new vertical position to particle
576
577	! ztold=ztra1(ipart)
578	zp=ztra1(ipart)
579
580	CALL redist_kf(lconvection,ldirect,delt,delx, & ! IN
581	dz1d,nzmax, nz,umf,uer,udr,dmf, & ! IN
582	der,ddr,rho1d, & ! IN
583	zf,zh, & ! IN
584	umfzf,dmfzf, & ! IN
585	zp) ! IN/OUT
586
587	if (zp .lt. 0.0) zp=-1.0*zp
588	if (zp .gt. height(nz)-0.5) &
589	zp = height(nz)-0.5
590	ztra1(ipart) = zp
591
592	! Calculate the gross fluxes across layer interfaces
593	!***************************************************
594
595	if (iflux.eq.1) then
596	itage=abs(itra1(ipart)-itramem(ipart))
597	do nage=1,nageclass
598	if (itage.lt.lage(nage)) goto 47
599	enddo
600	47 continue
601
602	if (nage.le.nageclass) &
603	call calcfluxes(nage,ipart,real(xtra1(ipart)), &
604	real(ytra1(ipart)),ztold)
605	endif
606
607	endif !(lconv .eqv. .true.)
608
609	enddo
610	60 continue
611	enddo
612	70 continue !inest - loop
613	! print*,'end of convmix'
614	!--------------------------------------------------------------------------
615	! write(,)'############################################'
616	! write(,)'TIME=',
617	! & itime
618	! write(,)'fraction of particles under convection',
619	! & sumconv/(sumall+0.001)
620	! write(,)'total number of particles',
621	! & sumall
622	! write(,)'number of particles under convection',
623	! & sumconv
624	! write(,)'############################################'
625
626	return
627	end subroutine convmix_kfeta

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format