Context Navigation

snuages3D.F @ 3094

Last change on this file since 3094 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 46.6 KB

Line
1	subroutine snuages(ngrid,tab1,tab2,tab3,tab4,tab5,
2	& x1,xnz,xnrad,ihor,fluxi,taused,precip)
3
4	!*
5	!
6	!1 c quantite de noyaux de la grille de rayon r a l'altitude z
7	!2 c quantite de glace 1 de la grille de rayon r a l'altitude z
8	!3 c quantite de glace 2 de la grille de rayon r a l'altitude z
9	!4 c quantite de glace 3 de la grille de rayon r a l'altitude z
10	!5 c quantite d'aerosols de la grille de rayon r a l'altitude z
11	c
12	c
13	!x1 timestep.
14
15	c Notes : taused n'est actuellement pas calculé !
16	c taused = dz/vitesse(z)
17	c
18	c
19	--------------------------------------------------------------
20	* *
21	* ENTRE 0 ET 1000 KILOMETRES *
22	* *
23	* la dimension fractale est en tableau, attention au *
24	* raccordement entre le regime moleculaire et le regime *
25	* fluide *
26	* *
27	* Modele microphysique: Cabane et al.,1992 / *
28	* Modele version fractale: Cabane et al.,1993 / *
29	* *
30	--------------------------------------------------------------
31	* VERSION DU 2 JUIN 1993 --- AUT 1994 --- 11/04/96
32	*
33	* changer: altitude de production z0=/taux de production ctot=
34	* : la charge/micron, ne
35	* : df(h),rf...
36	* raccordement aknc
37	*
38	* declaration des blocs communs
39	*------------------------------
40
41	use dimphy
42	IMPLICIT NONE
43	#include "dimensions.h"
44	#include "microtab.h"
45	#include "varmuphy.h"
46
47	common/x2ctps/li,lf,dt
48	common/x2con/c,c1,c2,c3,caer
49	common/cldpart/rad,mmu
50	common/x2frac/rfg,dfg
51	common/x2effets/ xsaison
52
53
54	* declaration des variables communes
55	* ----------------------------------
56
57	integer xnz,xnrad,ngrid
58	integer i,j,k,ihor,ibidon
59	integer ihormx,imx
60	integer li,lf
61	real xsaison
62	real dt
63	real c(nz,nrad,2)
64	real c1(nz,nrad,2),c2(nz,nrad,2)
65	real c3(nz,nrad,2)
66	real caer(nz,nrad,2)
67	real rad(nz,nrad), mmu(nz,nrad)
68	real rtemp(nz,nrad),mmutemp(nz,nrad)
69	real knu,knu2,w
70	real rfg(nz),dfg(nz,nrad)
71	real ddt,dtau,vitesse2
72	real vmax,vmin,rmax,rmin
73	real fluxi(ngrid,nz,3)
74	real taused(ngrid,nz)
75	real precip(ngrid,5)
76	real rmx
77
78
79	* variables internes
80	* ------------------
81
82	integer h,ti,itime,icoal
83	real tab1(nz,nrad),tab2(nz,nrad),
84	& tab3(nz,nrad),tab4(nz,nrad)
85	real tab5(nz,nrad)
86	real x1,xvolume,xmasse,xnoyaux
87	real knd1,knd2,knd3,knd4,knd5,knd6
88
89	real dice1,dice2,dice3,dice4
90
91	data itime/0/
92	data icoal/0/
93	save itime,icoal
94	save vmax,vmin,rmax,rmin
95
96
97	* controles
98	* ---------
99
100	if (nrad.ne.xnrad) stop 'nrad.ne.xnrad'
101	if (nz.ne.xnz) stop 'nz.ne.xnz'
102
103	do i=1,nz
104	do j=1,nrad
105	c(i,j,1)=tab1(i,j)
106	c(i,j,2)=0.0
107	c1(i,j,1)=tab2(i,j)
108	c1(i,j,2)=0.0
109	c2(i,j,1)=tab3(i,j)
110	c2(i,j,2)=0.0
111	c3(i,j,1)=tab4(i,j)
112	c3(i,j,2)=0.0
113	caer(i,j,1)=tab5(i,j)
114	caer(i,j,2)=0.0
115	enddo
116	enddo
117
118
119	* initialisation
120	* --------------
121
122	* Parametres physiques de Titan
123
124	if (itime.eq.0) then
125	call init2
126	print*,'init2 dans snuages.F'
127	itime=1
128	endif
129
130	aknc=2.92 !<--------Df=3
131
132	* Propriete des gouttes nuageuses/ calculee ici pour la sedimentation
133	* r(i,j)=r (cst pour la colonne). Cela evite les accumulations dans
134	* certains niveaux du a des differences de vitesse de sedimentation...
135
136	* ATTENTION, CES DEFINITIONS SONT INTERNES/ VOIR MUPHYS.F POUR DEFINITIONS
137	* EXTERNES (UTILISEES DANS OPTCV ET OPTCI, PAR EX.)
138
139	do i=1,nz
140	xnoyaux=0.
141	xvolume=0.
142	xmasse=0.
143	do j=1,nrad
144	dfg(i,j)=3.00
145	rfg(j)=6.6e-8 !<--- arbitraire pour df=3
146	xnoyaux=xnoyaux+c(i,j,1)
147	xvolume=xvolume+c1(i,j,1)+c2(i,j,1)+c3(i,j,1)+
148	& v_e(j)*c(i,j,1)
149	xmasse =xmasse+
150	& c1(i,j,1)*rhoi_ch4+
151	& c2(i,j,1)*rhoi_c2h6+
152	& c3(i,j,1)*rhoi_c2h2+
153	& v_e(j)c(i,j,1)rhol
154	enddo
155
156	do j=1,nrad
157
158	if (xnoyaux .le. 1.e-25 .or. xvolume.eq.0.) then ! utile ?
159	dfg(i,j)=3.00
160	rtemp(i,j) = 1.e-9
161	mmutemp(i,j)= rhol !on prend la masse vol des aerosols
162	else
163	c Mais pourquoi se compliquer la vie alors que de toute facon on
164	c reforce la masse volumique a son veritable calcul ^^
165	c
166	if ((c1(i,j,1)+c2(i,j,1)+c3(i,j,1))/xvolume.le.0.1) then ! glace / total
167	dfg(i,j)=3.00
168	rtemp(i,j) = ( (xvolume/xnoyaux)0.75/pi)*(1./3.)
169	c mmutemp(i,j) = 1000.
170	mmutemp(i,j)= xmasse/xvolume
171	if(rtemp(i,j).gt.3.e-4) rtemp(i,j)=3.e-4
172	else
173	dfg(i,j)=3.00
174	rtemp(i,j) = ( (xvolume/xnoyaux)0.75/pi)*(1./3.)
175	mmutemp(i,j)= xmasse/xvolume
176	if(rtemp(i,j).gt.3.e-4) rtemp(i,j)=3.e-4
177	endif
178	endif
179	enddo
180	enddo
181
182	do i=1,nz
183	do j=1,nrad
184	rad(i,j)=rtemp(i,1) ! mm valeur qqsoit j
185	mmu(i,j)=mmutemp(i,1)
186	enddo
187	enddo
188
189
190	* Coefficients de coagulation
191
192	if (icoal.eq.1) call calcoag2 ! 2 <-- 1
193	vmin=1.e+6
194	vmax=0.
195	rmin=1.e+6
196	rmax=0.
197	ihormx=1
198	imx=1
199
200
201	* choix interne du temps d iteration
202	* ----------------------------------
203
204	dt=x1
205
206	* iteration du modele sur le temps
207	* ---------------------------------
208
209	li=1
210	lf=2
211
212	* li=1 lf=2
213
214	call sedifn_fast(ihor,dice1,dice2,dice3,dice4) ! 1 <-- 1
215	! call sedifn ! 2 <-- 1
216
217	li=3-li ! li devient 2
218	lf=3-lf ! lf devient 1
219
220	if (icoal.eq.1) then
221	* li=2 lf=1
222	call coagul2(ihor) ! 1 <-- 2
223
224	li=3-li ! li devient 1
225	lf=3-lf ! lf devient 2
226
227	endif
228
229	do i=1,nz
230	vmin=vitesse2(i,1,1)
231	do j=1,nrad
232	fluxi(ihor,nz+1-i,1)=fluxi(ihor,nz+1-i,1)
233	& -vminc1(i,j,li)rhoi_ch4
234	fluxi(ihor,nz+1-i,2)=fluxi(ihor,nz+1-i,2)
235	& -vminc2(i,j,li)rhoi_c2h6
236	fluxi(ihor,nz+1-i,3)=fluxi(ihor,nz+1-i,3)
237	& -vminc3(i,j,li)rhoi_c2h2
238	enddo
239	enddo
240
241	c En theorie, les diceX sont NEGATIF (en sedimentant on ne fait que perdre de la glace)
242	c Les precipitations sont comptees positivement. (ET ON NE PREND QUE DES VALEURS POSITIVES !)
243
244	precip(ihor,1)=AMAX1(-dice1/rhoi_ch4,0.) ! m3/m2 = m :)
245	precip(ihor,2)=AMAX1(-dice2/rhoi_c2h6,0.)
246	precip(ihor,3)=AMAX1(-dice3/rhoi_c2h2,0.)
247	precip(ihor,4)=AMAX1(-dice4/rhol,0.)
248
249	do i=1,nz
250	do j=1,nrad
251	tab1(i,j)=c(i,j,li) ! li=1
252	tab2(i,j)=c1(i,j,li) ! li=1
253	tab3(i,j)=c2(i,j,li) ! li=1
254	tab4(i,j)=c3(i,j,li) ! li=1
255	enddo
256	enddo
257
258	return
259
260	end
261
262
263
264	*______________________________________________________________________
265
266	real function lambda2(j,indic)
267	*
268	------------------------------------------------------------------
269	* fonction calculant le libre parcours moyen des molecules *
270	* atmospheriques( rayon =ra) se trouvant dans la couche no j. *
271	* pour indic=0 ...... la particule se trouve a la frontiere entre*
272	* les couches j et j-1 *
273	* pour indic=1 ...... la particule se trouve au milieu de la *
274	* la couche j *
275	------------------------------------------------------------------
276	*
277	* declaration des blocs communs
278	*------------------------------
279	use dimphy
280	IMPLICIT NONE
281	#include "dimensions.h"
282	#include "microtab.h"
283	#include "varmuphy.h"
284
285	common/x2frac/rfg,dfg
286
287	* declaration des variables communes
288	* ----------------------------------
289
290	real rfg(nz),dfg(nz,nrad)
291
292	* declaration des variables internes
293	* ----------------------------------
294
295	integer indic,j
296	real pp,ra
297
298	ra=1.75e-10
299
300	* traitement
301	* ----------
302
303	if (indic.eq.0) then
304	pp=pb(j)
305	else
306	if (indic.ne.1) then
307	print*,'erreur argument fonction lambda'
308	stop
309	endif
310	pp=p(j)
311	endif
312
313	lambda2=kbzt(j)/(4sqrt(2.)pi(ra*2)pp)
314	end
315
316	*******************************************************************************
317
318	real function knu2(j,k,indic)
319	*
320	--------------------------------------------------------------
321	* fonction calculant le nombre de knudsen d'une particule *
322	* d'aerosol de rayon rad(k) se trouvant dans la couche no j *
323	* indic ...... idem function lambda *
324	--------------------------------------------------------------
325	*
326	* declaration des blocs communs
327	*------------------------------
328	use dimphy
329	IMPLICIT NONE
330	#include "dimensions.h"
331	#include "microtab.h"
332	#include "varmuphy.h"
333
334	common/cldpart/rad,mmu
335	common/x2frac/rfg,dfg
336
337
338	* declaration des variables communes
339	* ----------------------------------
340
341	real rad(nz,nrad),mmu(nz,nrad)
342	real rfg(nz),dfg(nz,nrad)
343
344
345	* declaration des variables internes
346	* ----------------------------------
347
348	integer indic,j,k
349	real lambda2,rfk
350
351	* traitement
352	* ----------
353
354	if (indic.ne.0 .and.indic.ne.1) then
355	print*,'erreur argument fonction knu'
356	stop
357	endif
358
359
360	rfk=(rad(j,k)*(3./dfg(j,k)))((rfg(k))**(1.-3./dfg(j,k)))
361	knu2=lambda2(j,indic)/rfk
362
363	end
364
365	*****************************************************************************
366
367	real function nud2(j,indic)
368	*
369	--------------------------------------------------------------
370	* fonction calculant la viscosite dynamique (en USI) de l air *
371	* d apres la formule de Sutherlant a l altitude j *
372	* indic ......... idem fonction lambda *
373	--------------------------------------------------------------
374	*
375	use dimphy
376	IMPLICIT NONE
377	#include "dimensions.h"
378	#include "microtab.h"
379	#include "varmuphy.h"
380	integer indic,j
381	real nud0,c,tt
382	real rfg(nz),dfg(nz,nrad)
383
384	common/x2frac/rfg,dfg
385
386	*
387	nud0=1.74e-5
388	c=109.
389
390	if(indic.ne.0.and.indic.ne.1) then
391	print*,'erreur argument fonction nud'
392	stop
393	endif
394
395	if (indic.eq.0) tt=tb(j)
396	if (indic.eq.1) tt=t(j)
397	nud2=nud0sqrt(tt/293)(1+c/293)/(1+c/tt)
398	end
399
400	****************************************************************************
401
402	real function vitesse2(j,k,indic)
403	*
404	-----------------------------------------------------------------
405	* fonction calculant la vitesse de chute d une particule de rayon*
406	* k se trouvant a l altitude j suivant la valeur du nombre de *
407	* Knudsen *
408	* indic ....... idem function lambda *
409	-----------------------------------------------------------------
410	*
411
412	* declaration des blocs communs
413	*------------------------------
414	use dimphy
415	IMPLICIT NONE
416	#include "dimensions.h"
417	#include "microtab.h"
418	#include "varmuphy.h"
419	common/cldpart/rad,mmu
420	common/x2frac/rfg,dfg
421	common/x2ctps/li,lf,dt
422
423	* declaration des variables communes
424	* ----------------------------------
425
426	real rad(nz,nrad),mmu(nz,nrad)
427	real rfg(nz),dfg(nz,nrad)
428
429
430	* declaration des variables internes
431	* ----------------------------------
432
433	integer indic,j,k,li,lf
434	real w,g,m,a0,zz,knu2,nud2,knud,tt,rhoh
435	real vlimite,akncx,rbis,rfk
436	real dt
437
438	* traitement
439	* ----------
440
441	if (indic.ne.0.and.indic.ne.1) then
442	print*,'erreur argument fonction vitesse'
443	stop
444	endif
445
446	if(indic.eq.0) then
447	zz=z(j)+dz(j)/2.
448	tt=tb(j)
449	rhoh=rhob(j)
450	endif
451	if(indic.eq.1) then
452	zz=z(j)
453	tt=t(j)
454	rhoh=rho(j)
455	endif
456
457	g=g0(rtit/(rtit+zz))*2
458	a0=0.74
459	m=(ach4(j)mch4+aar(j)mar+an2(j)*mn2)/nav
460	knud=knu2(j,k,indic)
461
462	c akncx=aknc
463	c if(df(k).gt.2.5) akncx=2.7
464
465	c if(knud.ge.akncx) then
466	c rbis=(rad(j,k)*(3.-6./dfg(j,k)))((rfg(k))**(-2.+6./dfg(j,k)))
467	c w=a0grbismmu(j,k)/(rhohsqrt(8kbztt/(pi*m)))
468	c endif
469
470	rfk=(rad(j,k)*(3./dfg(j,k)))((rfg(k))**(1.-3./dfg(j,k)))
471	w=2./9.rfk(dfg(j,k)-1.)rfg(k)*(3.-dfg(j,k))g*mmu(j,k)
472	& /nud2(j,indic)
473
474	w=w(1+1.2517knud+0.4knudexp(-1.1/knud))
475
476	w=w!*3. ! on tient compte de la largeur de distribution... a affiner
477	vitesse2=w
478
479
480	end
481	***********************************************************************
482	real function kd2(h)
483	*
484	--------------------------------------------------------------------
485	* cette fonction calcule le coefficient du terme de eddy diffusion *
486	* a l altitude j *
487	--------------------------------------------------------------------
488	*
489	use dimphy
490	IMPLICIT NONE
491	#include "dimensions.h"
492	#include "microtab.h"
493	#include "varmuphy.h"
494
495	common/x2frac/rfg,dfg
496
497	real zbx
498	real rfg(nz),dfg(nz,nrad)
499
500	integer h
501
502	zbx=z(h)+dz(h)/2.
503	if(zbx.le.42000.) then
504	c kd2=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
505	kd2=4.
506	else
507	kd2=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
508	kd2=0.0*kd2
509	endif
510	kd2=00.
511
512
513	return
514	end
515
516
517	*____________________________________________________________________________
518
519	subroutine init2
520	*
521	--------------------------------------------------------------------
522	* cette routine effectue : *
523	* 1) interpolation a partir des donnees initiales des *
524	* valeurs de p,t,rho,ach4,aar,an2 sur la grille *
525	* 2) initialisation des constantes (common/x2phys/) *
526	* 3) initialisation des variables temporelles (common *
527	* /temps/) *
528	* 4) definition des grilles en rayon et verticale *
529	* 5) initialisation de c(z,r,t) avec les donnees du *
530	* fichier unit=1 *
531	* *
532	* les donnees sont des valeurs caracterisques de l atmosphere de *
533	* TITAN ( voir Lelouch and co ) *
534	--------------------------------------------------------------------
535
536	* declaration des blocs communs
537	*------------------------------
538	use dimphy
539	IMPLICIT NONE
540	#include "dimensions.h"
541	#include "microtab.h"
542	#include "varmuphy.h"
543	common/x2ctps/li,lf,dt
544	common/cldpart/rad,mmu
545
546	* declaration des variables communes
547	* ----------------------------------
548
549	integer li,lf
550	real dt
551	real rad(nz,nrad), mmu(nz,nrad)
552
553
554	* declaration des variables internes
555	* ----------------------------------
556	integer nzd,i,ii
557	parameter (nzd=254)
558	integer limsup,liminf,j1,j2
559	real zd(nzd),ach4d(nzd),rap
560	real m
561
562
563	* initialisation des variables temporelles
564	* ----------------------------------------
565
566	li=1
567	lf=2
568
569
570	* interpolation de xch4,xar et xn2 sur la grille
571	* ----------------------------------------------
572
573	* donnees initiales (Lellouch et al,87)
574	* -------------------------------------
575
576	c print,'***** init'
577	do 1 i=1,168
578	zd(i)=(1000.-5(i-1))1000.
579	1 continue
580	do 2 i=1,78
581	zd(168+i)=(160.-2(i-1))1000.
582	2 continue
583	do 3 i=1,4
584	zd(246+i)=(5.-(i-1))*1000.
585	3 continue
586	do 4 i=1,4
587	zd(250+i)=(1.5-(i-1)0.5)1000.
588	4 continue
589
590	data (ach4d(i),i=1,168)/168*1.5e-2/
591	data (ach4d(i),i=169,254)/63*1.5e-2,1.6e-2,1.8e-2,1.8e-2,
592	& 1.9e-2,2.e-2,2.1e-2,2.3e-2,2.5e-2,2.8e-2,3.1e-2,3.6e-2,
593	& 4.1e-2,4.7e-2,5.7e-2,6.7e-2,7.5e-2,7*8.e-2/
594
595	liminf=0
596	limsup=0
597
598	* interpolation des taux de melange de ch4,ar,n2
599	*-----------------------------------------------
600
601	do 20 j1=1,nz
602	do 21 j2=1,nzd
603	if( zd(j2).le.z(j1)) goto 22
604	21 continue
605	22 liminf=j2
606	if (zd(liminf).eq.z(j1) )then
607	ach4(j1)=ach4d(liminf)
608	goto 20
609	endif
610	if (j2.ne.1) then
611	limsup=j2-1
612	else
613	limsup=j2
614	endif
615	if (limsup.eq.liminf) then
616	ach4(j1)=ach4(limsup)
617	else
618	ach4(j1)=ach4d(liminf)-(ach4d(limsup)-ach4d(liminf))/
619	s (zd(limsup)-zd(liminf))*(zd(liminf)-z(j1))
620	endif
621	20 continue
622
623	* rap= aar/an2 cst sur l altitude
624
625	rap=0.191
626	do 23 i=1,nz
627	an2(i)=(1.-ach4(i))/(1.+rap)
628	aar(i)=rap*an2(i)
629	23 continue
630
631	do 24 i=1,nz
632	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
633	rho(i)=p(i)m/(rgpt(i))
634	24 continue
635
636	do 34 i=1,nz
637	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
638	rhob(i)=pb(i)m/(rgptb(i))
639	c print*,pb(i),m,rgp,tb(i),rhob(i),rho(i)
640	34 continue
641
642	* fin d interpolation des taux de melange
643	*----------------------------------------
644
645	c print,'*** fin init'
646	540 continue
647	return
648
649	500 print*,'erreur lecture initialisation de c...erreur=',ii
650	stop
651
652	end
653
654	*__________________________________________________________________________
655
656	subroutine sedifn
657	*
658	------------------------------------------------------------------
659	* cette routine calcule l evolution de la fonction de distribution*
660	* c(z,r,t) pour les phenomenes de sedimentation et de diffusion *
661	------------------------------------------------------------------
662	*
663	*
664	* declaration des blocs communs
665	*------------------------------
666	use dimphy
667	IMPLICIT NONE
668	#include "dimensions.h"
669	#include "microtab.h"
670	#include "varmuphy.h"
671
672	common/x2ctps/li,lf,dt
673	common/x2con/ctmp,c1,c2,c3,caer
674	common/cldpart/rad,mmu
675	common/x2frac/rfg,dfg
676
677
678	* declaration des variables communes
679	* ----------------------------------
680
681	integer li,lf
682	real dt
683	real c(nz,nrad,2)
684	real c1(nz,nrad,2)
685	real c2(nz,nrad,2)
686	real c3(nz,nrad,2)
687	real caer(nz,nrad,2)
688	real ctmp(nz,nrad,2)
689	real rad(nz,nrad),mmu(nz,nrad)
690	real rfg(nz),dfg(nz,nrad)
691
692
693	* declaration des variables internes
694	* ----------------------------------
695
696	real w,w1,dzbX,dc
697	integer i,j,k,ii,nb
698	double precision sigma,theta,hc,l,rap,cmp,wp
699	double precision fs(nz+1),ft(nz+1)
700	real as(nz),bs(nz),cs(nz),ds(nz)
701	double precision asi(nztop:nz),bsi(nztop:nz),
702	& csi(nztop:nz)
703	double precision dsi(nztop:nz),xsol(nztop:nz)
704	real vitesse2,kd2
705
706	external dtridgl
707
708	* resolution
709	*------------
710
711	do 5 ii=1,4
712	do k=1,nrad
713	do j=nztop,nz
714	if( ii.eq.1 ) c(j,k,li)=ctmp(j,k,li)
715	if( ii.eq.2 ) c(j,k,li)=c1(j,k,li)
716	if( ii.eq.3 ) c(j,k,li)=c2(j,k,li)
717	if( ii.eq.4 ) c(j,k,li)=c3(j,k,li)
718	enddo
719	enddo
720
721	do 10 k=1,nrad
722	do 20 j=nztop,nz
723	if (j.eq.1) goto 20
724	* calcul de la vitesse corrigee
725
726	dzbX=(dz(j)+dz(j-1))/2.
727	w= -1*vitesse2(j,k,0)
728
729	if (kd2(j).ne.0.) then
730	theta=0.5(wdzbX/kd2(j)+log(rho(j-1)/rho(j)))
731	if (theta.ne.0) then
732	sigma=1./dtanh(theta)-1./theta
733	else
734	sigma=1.
735	endif
736	else
737	sigma=1.
738	endif
739
740	if(c(j,k,li).eq.0.) then
741	rap=10.
742	else
743	rap=c(j-1,k,li)/c(j,k,li)
744	if( rap.gt.10.) rap=10.
745	if( rap.lt.0.1) rap=0.1
746	endif
747
748	if (rap.gt.0.9 .and. rap.lt.1.1) then
749	w1=w
750	else
751	if(w.ne.0) then
752	hc=dzbX/dlog(rap)
753	l=dzbX/(wdt)(dexp(-w*dt/hc)-1.)/(1.-rap)
754	wp=w*1.d0
755	cmp=dlog(-wp)+abs(sigma)*dlog(l)
756	if (cmp.gt.38) then
757	goto 20
758	endif
759	w1=-dexp(cmp)
760	else
761	w1=0.
762	endif
763	endif
764
765	* calcul des flux aux interfaces
766
767
768	if (kd2(j).ne.0.) then
769	if (theta.ne.0.) then
770	ft(j)=(w1+log(rho(j-1)/rho(j))*kd2(j)/dzbX)/
771	& (dexp(2.*theta)-1.)
772	fs(j)=ft(j)dexp(2.theta)
773	else
774	ft(j)=kd2(j)/dzbX
775	fs(j)=kd2(j)/dzbX
776	endif
777	else
778	if (w1.lt.0.)then
779	ft(j)=-w1
780	fs(j)=0.
781	else
782	ft(j)=0.
783	fs(j)=w1
784	endif
785	endif
786
787	20 continue
788
789	* conditions aux limites pour les flux aux interfaces
790
791	fs(1)=0.
792	ft(1)=0.
793	fs(nz+1)=0.
794	ft(nz+1)=-w1
795
796	* calcul des coefficients de l equation discrete
797
798	do 21 j=nztop,nz
799	as(j)=-dz(j)/dt
800	bs(j)=-ft(j)
801	cs(j)=ft(j+1)+fs(j)-dz(j)/dt
802	ds(j)=-fs(j+1)
803	if ( cs(j).gt.0) goto 100
804	21 continue
805
806	* cas explicite (mu=0) : calcul de la fonction c(z,r,t+1)
807
808	do 22 j=nztop,nz-1
809
810	if (j.eq.nztop) then
811	dc=(cs(nztop)*c(nztop,k,li)+ds(nztop)
812	& *c(nztop+1,k,li))/as(nztop)
813	c(nztop,k,lf)=dc
814	goto 22
815	endif
816
817	dc=(bs(j)c(j-1,k,li)+cs(j)c(j,k,li)+ds(j)*c(j+1,k,li))
818	s /as(j)
819	c(j,k,lf)=dc
820
821	22 continue
822
823	dc=(bs(nz)c(nz-1,k,li)+cs(nz)c(nz,k,li))/as(nz)
824	c(nz,k,lf)=dc
825
826	if (nztop.ne.1) then
827	do 32 j=1,nztop-1
828	c(j,k,lf)=c(j,k,li)
829	32 continue
830	endif
831
832	goto 10
833
834	100 continue
835
836	* cas implicite (mu=1) : calcul de la fonction c(z,r,t+1)
837
838	do 101 j=nztop,nz
839	asi(j)=ft(j)
840	bsi(j)=-(ft(j+1)+fs(j)+dz(j)/dt)
841	csi(j)=fs(j+1)
842	dsi(j)=-dz(j)/dt*c(j,k,li)
843	101 continue
844
845	* inversion de la matrice tridiagonale
846
847	nb=nz-nztop+1
848
849	call dtridgl(nb,asi,bsi,csi,dsi,xsol)
850
851	do 102 j=nztop,nz
852	c(j,k,lf)=xsol(j)
853	102 continue
854
855	if (nztop.ne.1) then
856	do 110 j=1,nztop-1
857	c(j,k,lf)=c(j,k,li)
858	110 continue
859	endif
860
861
862
863	10 continue
864
865	do k=1,nrad
866	do j=nztop,nz
867	if( ii.eq.1 ) ctmp(j,k,lf)=c(j,k,lf)
868	if( ii.eq.2 ) c1(j,k,lf) =c(j,k,lf)
869	if( ii.eq.3 ) c2(j,k,lf) =c(j,k,lf)
870	if( ii.eq.4 ) c3(j,k,lf) =c(j,k,lf)
871	enddo
872	enddo
873
874	5 continue
875
876	return
877
878	end
879
880	*__________________________________________________________________________
881
882	subroutine sedifn_fast(ihor,dice1,dice2,dice3,dice4)
883	*
884	------------------------------------------------------------------
885	* cette routine calcule l evolution de la fonction de distribution *
886	* c(z,r,t) pour les phenomenes de sedimentation {pas de diffusion} *
887	* dice1 = delta glace CH4 *
888	* dice2 = delta glace C2H6 *
889	* dice3 = delta glace C2H2 *
890	* dice4 = delta Volume noyaux *
891	------------------------------------------------------------------
892	*
893	* declaration des blocs communs
894	*------------------------------
895	use dimphy
896	IMPLICIT NONE
897	#include "dimensions.h"
898	#include "microtab.h"
899	#include "varmuphy.h"
900
901	common/x2ctps/li,lf,dt
902	common/x2con/c,c1,c2,c3,caer
903	common/cldpart/rad,mmu
904	common/x2frac/rfg,dfg
905
906
907	* declaration des variables communes
908	* ----------------------------------
909
910	integer li,lf,i,j,k,ihor
911	integer jinf,jsup,jj,iiter
912	real dt
913	real c(nz,nrad,2)
914	real c1(nz,nrad,2)
915	real c2(nz,nrad,2)
916	real c3(nz,nrad,2)
917	real caer(nz,nrad,2)
918	real ci(nz,nrad,2)
919	real ci1(nz,nrad,2)
920	real ci2(nz,nrad,2)
921	real ci3(nz,nrad,2)
922	real rad(nz,nrad),mmu(nz,nrad)
923	real rfg(nz),dfg(nz,nrad)
924	real puit(nz)
925	c ------ echange est cree sur la taille maxi mais n'est utilisee
926	c que sur la dim geree par le proc (klon ou jjm+1)
927	integer ngrid
928	parameter (ngrid=(jjm-1)*iim+2) ! = taille maximum
929	real echange(nz,nz,ngrid)
930	c pas genial mais vu que c est tres local, pas de soucis a priori en parallele.
931	real bilan1,bilan2,bilan3,bilan4,bilan5
932	real bilan11,bilan12,bilan13,bilan14,bilan15
933	real compte,compte2,xepl
934	real dice1,dice2,dice3,dice4
935
936
937
938	* declaration des variables internes
939	* ----------------------------------
940
941	real vitesse2,kd2
942	real w,ws,wi,zs,zi,alpha,v0,deltazs,deltazi
943	real zni,znip1,xcnt,ichx,arg1,arg2,xft,xf
944	real argexp,explim
945
946	save echange
947
948
949	* resolution
950	*------------
951
952
953	bilan1=0.
954	bilan2=0.
955	bilan3=0.
956	bilan4=0.
957	bilan5=0.
958	do k=1,nrad
959	do j=nztop,nz
960	ci(j,k,li)= c(j,k,li)*dzb(j) ! li
961	ci1(j,k,li)=c1(j,k,li)*dzb(j)
962	ci2(j,k,li)=c2(j,k,li)*dzb(j)
963	ci3(j,k,li)=c3(j,k,li)*dzb(j)
964	bilan5=bilan5+ci(j,k,li)
965	bilan1=bilan1+ci1(j,k,li)
966	bilan2=bilan2+ci2(j,k,li)
967	bilan3=bilan3+ci3(j,k,li)
968	bilan4=bilan4+
969	& ci(j,k,li)4./3.pirf(k)3.vrat_e**(k-imono)
970
971	ci(j,k,lf)= 0. ! lf
972	ci1(j,k,lf)=0.
973	ci2(j,k,lf)=0.
974	ci3(j,k,lf)=0.
975	enddo
976	enddo
977
978	* calcul de la matrice d echange
979	*----------------------------------------------------------------
980
981	do j=nztop,nz
982	do i=nztop,nz
983	echange(i,j,ihor)=0.
984	enddo
985	enddo
986
987	do 20 i=nztop,nz
988	puit(i)=0.
989	do 30 k=1,1
990	xcnt=0.
991	ICHX=1 ! extrapolation 0: lineaire 1: exponentielle
992	IF(ICHX.eq.0 .or. ICHX.eq.2) THEN
993	ws=vitesse2(i,k,0)
994	if(i.lt.nz) wi=vitesse2(i+1,k,0)
995	if(i.eq.nz) wi=vitesse2(i ,k,1)
996	w=(ws+wi)/2.
997	zni =zb(i)-w*dt
998	znip1=zb(i)-dzb(i)-w*dt
999	ENDIF
1000
1001	explim=30.
1002
1003	IF(ICHX.eq.1 .or. ICHX.eq.2) THEN
1004	ws=vitesse2(i,k,0)
1005	zs=zb(i)
1006	wi=vitesse2(i,k,1)
1007	zi=z(i)
1008
1009	c if(wi.eq.ws) wi=ws/1.001 ! sinon ca plante !
1010	if(abs(wi-ws)/wi .le. 1.e-3) wi=ws/1.001 ! sinon ca plante !
1011
1012	if(wi.ne.0.) alpha= alog(ws/wi) /(zs-zi)
1013	argexp=alpha*zs
1014	if(argexp.lt.-explim) argexp=-explim
1015	if(argexp.gt. explim) argexp=+explim
1016	v0 = ws/exp(argexp)
1017
1018	argexp=alpha*zb(i)
1019	if(argexp.lt.-explim) argexp=-explim
1020	if(argexp.gt. explim) argexp=+explim
1021	arg1=1.+v0alphaexp(argexp)*dt
1022
1023	argexp=alpha*(zb(i)-dzb(i))
1024	if(argexp.lt.-explim) argexp=-explim
1025	if(argexp.gt. explim) argexp=+explim
1026	arg2=1.+v0alphaexp(argexp)*dt
1027
1028	iiter=0
1029	101 continue
1030
1031	if (iiter.le.25) then
1032	if(arg1.le.0. .or. arg2.le.0.) then
1033	print*,ihor,i,iiter, 'ajustement vitesse',arg1,arg2
1034	print*,ws,wi, ' w1 w2 anc valeurs'
1035	print*,alpha, ' alpha anc valeurs'
1036	print*,rad(i,k), 'r(j,k)'
1037	print*,mmu(i,k), ' mmu(j,k)'
1038	print*,t(i), ' t(j,k)'
1039	print*,arg1,arg2, ' arg1 arg2 anc valeurs'
1040	arg2=1+(arg2-1.)/2.
1041	arg1=1+(arg1-1.)/2.
1042	iiter=iiter+1
1043	print*,arg1,arg2, ' arg1 arg2 nle valeurs'
1044	goto 101
1045	endif
1046	else
1047	stop
1048	endif
1049
1050	if(arg1.lt.0.) print*,'arg1:',arg1
1051	if(arg2.lt.0.) print*,'arg2:',arg2
1052
1053	deltazs=-alog(arg1)/alpha
1054	deltazi=-alog(arg2)/alpha
1055
1056	zni =zb(i)+deltazs
1057	znip1=zb(i)-dzb(i)+deltazi
1058
1059	ENDIF
1060
1061	if(zni.ne.znip1) xft=zni/(zni-znip1)
1062	if(zni.eq.znip1 .and. zni.le.0.) xft=0.
1063	if(zni.eq.znip1 .and. zni.gt.0.) then
1064	xft=0.
1065	print*,'zni..znip1', zni,znip1
1066	endif
1067
1068	* Si des aerosols touchent le sol (zni < 0 ) alors on fixe
1069	* le niveau a 0, et on elimine les aerosols correspondants
1070	*-----------------------------------------------------------
1071
1072	if(znip1 .lt. 0.) znip1=0.
1073	if(zni .lt. 0.) zni =0.
1074
1075	if(zni.le.0. .and. znip1 .le. 0.) then
1076	c print*,'voie 1 / disparition complete'
1077	xft=0.
1078	xf=1.
1079	xcnt=xcnt+xf
1080	puit(i)=puit(i)+xf
1081	endif
1082
1083	if(zni.gt.0. .and. znip1 .le. 0.) then
1084	c print*,'voie 2 / disparitipon partielle'
1085	xf=(1-xft)
1086	xcnt=xcnt+xf
1087	puit(i)=puit(i)+xf
1088	endif
1089
1090	if(zni.gt.0. .and. znip1 .gt. 0.) then
1091	c print*,'voie 3 / pas de disparition'
1092	xft=1.
1093	xf=(1-xft)
1094	xcnt=xcnt+xf
1095	puit(i)=puit(i)+xf
1096	endif
1097
1098	jsup=nz+1
1099	jinf=nz+1
1100	do j=nztop,nz
1101	if(zni.le.zb(j) .and. zni.ge.zb(j)-dzb(j)) jsup=j
1102	if(znip1.le.zb(j).and. znip1.ge.zb(j)-dzb(j)) jinf=j
1103	enddo
1104	if(zni .ge. 0. .and. zni .lt. 1.e-3) jsup=nz
1105	if(znip1 .ge. 0. .and. znip1 .lt. 1.e-3) jinf=nz
1106
1107
1108	* Volume inclu dans un seul niveau
1109	*----------------------------------
1110
1111	if (jsup .eq. jinf. and. jsup.ge. nz+1) then
1112	xcnt=xcnt+1.
1113	print*,'cas impossible'
1114	print*,'alpha= ',alpha
1115	print*,'ws wi ',ws,wi
1116	print*,'deltazs deltazi ',deltazs,deltazi
1117	print*,' r(i,k) mmu(i,k)', rad(i,k),mmu(i,k)
1118	print*,' t(j,k)',t(i)
1119	print*,zni,znip1,jsup,jinf
1120	print*,'STOP'
1121	STOP
1122	endif
1123
1124	if (jsup .eq. jinf. and. jsup.le. nz) then
1125	xf=1.
1126	xcnt=xcnt+xft*xf
1127	if(jinf.le.nz) then
1128	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1129	endif
1130	endif
1131
1132	* Volume a cheval sur 2 niveaux
1133	*--------------------------------
1134
1135	if (jinf .eq. jsup+1) then
1136	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
1137	xcnt=xcnt+xf*xft
1138	if(jsup.le.nz) then
1139	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
1140	endif
1141	xf=(zb(jinf)-znip1)/(zni-znip1)
1142	xcnt=xcnt+xf*xft
1143	if(jinf.le.nz) then
1144	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1145	endif
1146	endif
1147
1148	* Volume a cheval sur 3 ou plus de niveaux
1149	*------------------------------------------
1150
1151	if (jinf .gt. jsup+1) then
1152
1153	c print*,' voie C / dans N cases'
1154	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
1155	xcnt=xcnt+xf*xft
1156	if(jsup.le.nz) then
1157	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
1158	endif
1159
1160	xf=(zb(jinf)-znip1)/(zni-znip1)
1161	xcnt=xcnt+xf*xft
1162	if(jinf.le.nz) then
1163	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1164	endif
1165
1166	do jj=jsup+1,jinf-1
1167	xf=(zb(jj)-zb(jj+1))/(zni-znip1)
1168	xcnt=xcnt+xf*xft
1169	if(jj.le.nz) then
1170	echange(jj,i,ihor)=echange(jj,i,ihor)+xft*xf
1171	endif
1172	enddo
1173
1174	endif
1175
1176
1177	* et sur les rayons...
1178	*---------------------
1179	30 continue
1180
1181	* fin de la grande boucle sur les altitudes...
1182	*----------------------------------------------
1183	20 continue
1184
1185	* Calcul etat final Cfinal = Echange*initial
1186	*----------------------------------------------
1187
1188	compte=0.
1189	compte2=0.
1190	do j=1,nz
1191	xepl=0.
1192	do jj=1,nz
1193	xepl=xepl+echange(jj,j,ihor)
1194	compte=compte+echange(jj,j,ihor)
1195	compte2=compte2+echange(jj,j,ihor)
1196	enddo
1197	compte2=compte2+puit(j)
1198	enddo
1199
1200	if(abs(compte2-nz) .gt. 1.e-4)
1201	& print*,'Matrice calculee#',ihor,'tx d expl (/55):',
1202	& compte,compte2
1203
1204
1205	* Fin du calcul de la matrice d*echange
1206	*----------------------------------------------
1207
1208	do j=nztop,nz
1209	do k=1,nrad
1210
1211	do jj=nztop,nz
1212	ci(j,k,lf)=ci(j,k,lf)+
1213	& echange(j,jj,ihor)*ci(jj,k,li)
1214	ci1(j,k,lf)=ci1(j,k,lf)+
1215	& echange(j,jj,ihor)*ci1(jj,k,li)
1216	ci2(j,k,lf)=ci2(j,k,lf)+
1217	& echange(j,jj,ihor)*ci2(jj,k,li)
1218	ci3(j,k,lf)=ci3(j,k,lf)+
1219	& echange(j,jj,ihor)*ci3(jj,k,li)
1220	enddo
1221	enddo
1222	enddo
1223
1224
1225	* Controles et affichage Bilan
1226	*----------------------------------------------
1227
1228	bilan11=0.
1229	bilan12=0.
1230	bilan13=0.
1231	bilan14=0.
1232	bilan15=0.
1233	do k=1,nrad
1234	do j=nztop,nz
1235	c(j,k,lf) =ci(j,k,lf)/dzb(j)
1236	c1(j,k,lf)=ci1(j,k,lf)/dzb(j)
1237	c2(j,k,lf)=ci2(j,k,lf)/dzb(j)
1238	c3(j,k,lf)=ci3(j,k,lf)/dzb(j)
1239	bilan15=bilan15+ ci(j,k,lf)
1240	bilan11=bilan11+ci1(j,k,lf)
1241	bilan12=bilan12+ci2(j,k,lf)
1242	bilan13=bilan13+ci3(j,k,lf)
1243	bilan14=bilan14+
1244	& ci(j,k,lf)4./3.pirf(k)3.vrat_e**(k-imono)
1245	enddo
1246	enddo
1247
1248	c print*,'sedifn_fast Bilans:'
1249	c & ,bilan11,bilan12,bilan13
1250	c print*,'Bilan1:',bilan1,bilan11
1251	c print*,'Bilan2:',bilan2,bilan12
1252	c print*,'Bilan3:',bilan3,bilan13
1253	c print*,'Bilan5:',bilan5,bilan15
1254
1255	dice1=0.
1256	dice2=0.
1257	dice3=0.
1258	dice4=0.
1259	dice1=(bilan11-bilan1)rhoi_ch4 !glace 1 m^3.m^-2 kg.m^-3 pourquoi ????
1260	dice2=(bilan12-bilan2)*rhoi_c2h6 !glace 2
1261	dice3=(bilan13-bilan3)*rhoi_c2h2 !glace 3
1262	dice4=(bilan14-bilan4)*rhol !noyaux
1263
1264	return
1265	end
1266
1267
1268	subroutine coagul2(ihor)
1269
1270	*********************************************************
1271	* ce programme calcule la nouvelle concentration dans *
1272	* le a ieme intervalle de rayon, a l'altitude h, a *
1273	* l'instant t+dt *
1274	*********************************************************
1275
1276	use dimphy
1277	IMPLICIT NONE
1278	#include "dimensions.h"
1279	#include "microtab.h"
1280	#include "varmuphy.h"
1281
1282	* declaration des blocs communs
1283	*------------------------------
1284
1285	common/x2ctps/li,lf,dt
1286	common/x2con/c,c1,c2,c3,caer
1287	common/cldpart/rad,mmu
1288
1289	* declaration des variables
1290	* --------------------------
1291
1292	integer li,lf
1293	real dt
1294	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1295	& c3(nz,nrad,2)
1296	real caer(nz,nrad,2)
1297	real rad(nz,nrad),mmu(nz,nrad)
1298
1299	* declaration des variables propres au ss-programme
1300	* -------------------------------------------------
1301
1302	integer h,a,ihor,i
1303	real pr,pe,eta1,eta2
1304	real sum11,sum12,sum13,sum21,sum22,sum23
1305	real rfb,rfx,rpr
1306
1307	* traitement
1308	* ----------
1309
1310
1311	sum11=0.
1312	sum12=0.
1313	sum13=0.
1314	sum21=0.
1315	sum22=0.
1316	sum23=0.
1317
1318	do 10 h=nztop,nz
1319	do 11 a=1,nrad ! boucle aerosol secs
1320	call pertpro2(ihor,h,a,pe,pr)
1321
1322	if(1+dtpe .gt. 1.1) print,a,1+dt*pe,' scav'
1323	caer(h,a,lf)=(caer(h,a,li)+prdt)/(1+dtpe)
1324	c(h,a,lf)=c(h,a,li)
1325	c1(h,a,lf)=c1(h,a,li)
1326	c2(h,a,lf)=c2(h,a,li)
1327	c
1328	c eta1=0.
1329	c eta2=0.
1330	c if(c(h,a,li) .ne. 0.) eta1=c1(h,a,li)/c(h,a,li)
1331	c if(c(h,a,li) .ne. 0.) eta2=c2(h,a,li)/c(h,a,li)
1332	c c(h,a,lf) =( c(h,a,li)+prdt)/(1+dtpe)
1333	c c1(h,a,lf)=(c1(h,a,li)+eta1prdt)/(1+dt*pe)
1334	c c2(h,a,lf)=(c2(h,a,li)+eta2prdt)/(1+dt*pe)
1335	c
1336	c
1337	c sum11=sum11+c(h,a,li)
1338	c sum12=sum12+c1(h,a,li)
1339	c sum13=sum13+c2(h,a,li)
1340	c
1341	c sum21=sum21+ c(h,a,lf)
1342	c sum22=sum22+c1(h,a,lf)
1343	c sum23=sum23+c2(h,a,lf)
1344	c
1345	c
1346	11 continue
1347	10 continue
1348
1349
1350	if (nztop.ne.1) then
1351	do 12 h=1,nztop-1
1352	do 12 a=1,nrad
1353	c(h,a,lf)=c(h,a,li)
1354	c1(h,a,lf)=c1(h,a,li)
1355	c2(h,a,lf)=c2(h,a,li)
1356	caer(h,a,lf)=caer(h,a,li)
1357	12 continue
1358	endif
1359
1360	return
1361	end
1362
1363
1364	*__________________________________________________________________________
1365
1366	subroutine calcoag2
1367
1368	***************************************************************
1369	* *
1370	* Ce programme calcule les coefficients de collection d'une *
1371	* particule de rayon x avec une particule de rayon b a une *
1372	* altitude donnee h *
1373	***************************************************************
1374
1375	* declaration des blocs communs
1376	*------------------------------
1377	use dimphy
1378	IMPLICIT NONE
1379	#include "dimensions.h"
1380	#include "microtab.h"
1381	#include "varmuphy.h"
1382
1383	common/x2ctps/li,lf,dt
1384	common/x2con/c,c1,c2,c3,caer
1385	common/cldpart/rad,mmu
1386	common/x2coag/k
1387	common/x2frac/rfg,dfg
1388
1389	* declaration des variables
1390	* --------------------------
1391
1392	integer li,lf
1393	real dt
1394	real knu2,nud2,k(nz,nrad,nrad)
1395	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1396	& c3(nz,nrad,2)
1397	real caer(nz,nrad,2)
1398	real rad(nz,nrad),mmu(nz,nrad)
1399	real rfg(nz),dfg(nz,nrad)
1400
1401
1402	* declaration des variables propres au ss-programme
1403	* -------------------------------------------------
1404
1405	integer h,b,x,ihor,i
1406	real nua,lambb,lambx,knb,knx,alphab,alphax,d,e,f,kcg
1407	real db,dx,rm,dm,deltab,deltax,del,g,beta,gx,gb
1408	real rfb,rfx,rpr
1409	real*8 ne,qe,epso
1410	real*8 corelec,yy
1411
1412	real kco,vx,vb,vitesse2,sto,ee,a,dd,bb,p0,t0,l0,ccol
1413	real st(37),ef(37)
1414	real vitesse,knu
1415	external vitesse,knu
1416
1417
1418	* initialisation
1419	* --------------
1420
1421
1422
1423	* -nombres de STOCKES
1424
1425	data(st(i),i=1,37)/1.35,1.5,1.65,1.85,2.05,2.25,2.5,2.8,3.1,
1426	s 3.35,3.6,3.95,4.3,4.7,5.05,5.45,5.9,6.4,7.,7.6,8.3,9.05,9.9,
1427	s 10.9,11.1,13.5,15.3,17.25,20.5,24.5,30.4,39.3,48,57,86.,
1428	s 187.,600./
1429
1430	* -coef. d'efficacite de collection
1431
1432	ef(1)=3.75
1433	ef(2)=8.75
1434	do 11 i=3,37
1435	ef(i)=ef(i-1)+2.5
1436	11 continue
1437
1438	do 2 i=1,37
1439	ef(i)=ef(i)*1e-2
1440	2 continue
1441
1442	qe=1.6e-19
1443	ne=-30e+6
1444	epso=1e-9/(36*pi)
1445
1446	d=1.257
1447	e=0.4
1448	f=-1.1
1449
1450	* iteration sur z
1451
1452	do 1 h=1,nz
1453
1454	nua=nud2(h,1)
1455
1456	* iteration sur les rayons
1457
1458	do 1 b=1,nrad ! boucle aerosols secs : indice ''
1459
1460	knb=knu(h,b,1) ! knu et vitesse se trouvent dans brume.F
1461	vb=vitesse(h,b,1) ! ils concernent les aerosols
1462
1463	do 1 x=1,nrad !boucles gouttes : indice '2'
1464
1465	knx=knu2(h,x,1)
1466	vx=vitesse2(h,x,1)
1467
1468	COAGULATION **************************************************
1469	--------------**************************************************
1470	* calcul du terme correcteur 'slip-flow'
1471
1472	alphab=d+e*exp(f/knb)
1473	alphax=d+e*exp(f/knx)
1474
1475	* calcul du coefficient de diffusion
1476
1477
1478	rfb=(r_e(b)*(3./df(b))) ((rf(b)) **(1.-3./df(b)))
1479	rfx=(rad(h,x)**
1480	& (3./dfg(h,x)))((rfg(x))*(1.-3./dfg(h,x)))
1481
1482	db=kbzt(h)(1+alphabknb)/(6pinuarfb)
1483	dx=kbzt(h)(1+alphaxknx)/(6pinuarfx)
1484
1485	* calcul du coefficient de coagulation
1486
1487	rpr=rfb+rfx
1488	kcg=4pirpr*(db+dx)
1489
1490	* calcul de la vitesse thermique
1491
1492	gb=sqrt(6kbzt(h)/(rholpi2r_e(b)**3))
1493	gx=sqrt(6kbzt(h)/(rholpi2rad(h,x)**3))
1494
1495	* calcul du libre parcours apparent des aerosols
1496
1497	lambb=8db/(pigb)
1498	lambx=8dx/(pigx)
1499
1500	*calcul du terme correcteur beta
1501
1502	rm=rpr/2.
1503	dm=(dx+db)/2.
1504	g=sqrt(gx2+gb2)
1505	deltab=(((2rfb+lambb)3-(4rfb2+lambb2)**1.5)
1506	s /(6rfblambb)-2rfb)sqrt(2.)
1507	deltax=(((2rfx+lambx)3-(4rfx2+lambx2)**1.5)
1508	s /(6rfxlambx)-2rfx)sqrt(2.)
1509	del=sqrt(deltab2+deltax2)
1510	beta=1/((rm/(rm+del/2))+(4dm/(grm)))
1511
1512	* calcul du coefficient de coagulation corrige
1513
1514	kcg=kcg*beta
1515
1516	COALESCENCE ************************************************
1517	-------------************************************************
1518
1519	kco=0.
1520
1521	if ( b.eq. x) continue ! goto 9
1522
1523
1524	* calcul du nombre de Stockes de la petite particule
1525
1526	sto=2rholrfx*2abs(vx-vb)/(9nuarfb)
1527
1528	* calcul du coef. de Cunningham-Millikan
1529
1530	a=1.246
1531	bb=0.42
1532	dd=0.87
1533	l0=0.653e-7
1534	p0=101325.
1535	t0=288.
1536
1537	ee=1+
1538	& (l0t(h)p0(a+bbexp(-ddrfxt0p(h)/(l0t(h)*p0))))
1539	s /(rfxt0p(h))
1540
1541	* calcul du nombre de Stockes corrige
1542
1543	sto=sto*ee
1544
1545	if (sto .le. 1.2) goto 9
1546
1547	if (sto .ge. 600.) then
1548	ccol=1.
1549	goto 8
1550	endif
1551
1552	* recherche du coefficient de collection
1553
1554	do 3 i=1,37
1555	if (sto .gt. st(i)) then
1556	goto 3
1557	endif
1558	if (sto .eq. st(i)) then
1559	ccol=ef(i+1)
1560	else
1561	ccol=ef(i)
1562	endif
1563	goto 8
1564	3 continue
1565
1566	* calcul du coefficient de coalescence
1567
1568	8 kco=pi(rfb+rfx)2ccol*abs(vb-vx)
1569
1570	9 continue
1571
1572	CORRECTION ELECTRICITE *****************************
1573	------------------------****************************
1574
1575	c yy=1.d0ne2r(x)r(b)qe**2
1576	c & /(1.d0kbzt(h)(r(b)+r(x))4piepso)
1577	c corelec=0.
1578	c if (yy.lt.50.) corelec=yy/(exp(yy)-1.)
1579
1580	corelec=1.
1581
1582	c b=aerosol
1583	c x=gouttes
1584
1585	k(h,b,x)=(kcg+kco)*corelec
1586
1587	c
1588	c ATTENTION, IL N'Y A PLUS DE SYMETRIE...
1589	c k(ihor,h,x,b)=k(ihor,h,b,x)
1590	c
1591	c
1592
1593
1594	1 continue
1595	return
1596	end
1597
1598	*______________________________________________________________________
1599
1600	subroutine pertpro2(ihor,h,a,l_,pr_)
1601
1602	*****************************************************************************
1603	* *
1604	* ce programme permet le calcul du terme de production (pr) et de perte (l)*
1605	* pour le phenomene de coagulation *
1606	* dans le a ieme intervalle de rayon a une altitude h *
1607	****************************************************************************
1608
1609	* declaration des blocs communs
1610	*------------------------------
1611	use dimphy
1612	IMPLICIT NONE
1613	#include "dimensions.h"
1614	#include "microtab.h"
1615	#include "varmuphy.h"
1616
1617	common/x2ctps/li,lf,dt
1618	common/x2con/c,c1,c2,c3,caer
1619	common/cldpart/rad,mmu
1620	common/x2coag/k
1621
1622
1623	* declaration des variables
1624	* --------------------------
1625
1626	integer li,lf
1627	real dt
1628	real dr(nrad),dv(nrad)
1629	real k(nz,nrad,nrad)
1630	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1631	& c3(nz,nrad,2)
1632	real caer(nz,nrad,2)
1633	real rad(nz,nrad),mmu(nz,nrad)
1634
1635
1636	* declaration des variables propres au ss-programme
1637	* -------------------------------------------------
1638
1639	integer h,b,a,x,ihor,i
1640	real*8 pr,ss,s,l
1641	real pr_,l_,vol,del
1642
1643	* traitement
1644	* -----------
1645
1646
1647
1648	* production
1649	*+++++++++++++
1650	s=0.d0
1651	ss=0.d0
1652	pr=0.d0
1653
1654	c Pas de production d'aerosols par scavenging !!!
1655
1656	pr=0.d0
1657
1658	* perte
1659	*- - - - -
1660
1661	l=0.d0
1662
1663
1664	do 10 x=1,nrad ! boucle sur les gouttes
1665	l=l+k(h,a,x)*c(h,x,li)
1666	if(l.ne.0.d0) then
1667	print*,a,x,k(h,a,x),c(h,x,li),l,' : detail coal'
1668	endif
1669	10 continue
1670
1671
1672	#ifdef CRAY
1673	l_=l
1674	pr_=pr
1675	#else
1676	l_=sngl(l)
1677	pr_=sngl(pr)
1678	#endif
1679
1680	return
1681
1682	end
1683

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

Download in other formats:

Original Format