Context Navigation

snuages3D.F @ 201

Last change on this file since 201 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
Quelques ajouts dans la doc.

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

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

Download in other formats:

Original Format