Context Navigation

snuages3D.F @ 881

Last change on this file since 881 was 814, checked in by slebonnois, 12 years ago
SL: corrections de bugs pour Titan, + suppression de la dependance circulaire dans mod_grid_phy_lmdz.F90, + petite inversion sur la ligne 101 de create_make_gcm pour que disvert compile du premier coup.
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 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 decrit sur ngrid=klon mais peut etre utilisee
926	c uniquement sur jjm+1
927	integer ngrid
928	parameter (ngrid=(jjm-1)*iim+2) ! = klon
929	real echange(nz,nz,ngrid)
930	real bilan1,bilan2,bilan3,bilan4,bilan5
931	real bilan11,bilan12,bilan13,bilan14,bilan15
932	real compte,compte2,xepl
933	real dice1,dice2,dice3,dice4
934
935
936
937	* declaration des variables internes
938	* ----------------------------------
939
940	real vitesse2,kd2
941	real w,ws,wi,zs,zi,alpha,v0,deltazs,deltazi
942	real zni,znip1,xcnt,ichx,arg1,arg2,xft,xf
943	real argexp,explim
944
945	save echange
946
947
948	* resolution
949	*------------
950
951
952	bilan1=0.
953	bilan2=0.
954	bilan3=0.
955	bilan4=0.
956	bilan5=0.
957	do k=1,nrad
958	do j=nztop,nz
959	ci(j,k,li)= c(j,k,li)*dzb(j) ! li
960	ci1(j,k,li)=c1(j,k,li)*dzb(j)
961	ci2(j,k,li)=c2(j,k,li)*dzb(j)
962	ci3(j,k,li)=c3(j,k,li)*dzb(j)
963	bilan5=bilan5+ci(j,k,li)
964	bilan1=bilan1+ci1(j,k,li)
965	bilan2=bilan2+ci2(j,k,li)
966	bilan3=bilan3+ci3(j,k,li)
967	bilan4=bilan4+
968	& ci(j,k,li)4./3.pirf(k)3.vrat_e**(k-imono)
969
970	ci(j,k,lf)= 0. ! lf
971	ci1(j,k,lf)=0.
972	ci2(j,k,lf)=0.
973	ci3(j,k,lf)=0.
974	enddo
975	enddo
976
977	* calcul de la matrice d echange
978	*----------------------------------------------------------------
979
980	do j=nztop,nz
981	do i=nztop,nz
982	echange(i,j,ihor)=0.
983	enddo
984	enddo
985
986	do 20 i=nztop,nz
987	puit(i)=0.
988	do 30 k=1,1
989	xcnt=0.
990	ICHX=1 ! extrapolation 0: lineaire 1: exponentielle
991	IF(ICHX.eq.0 .or. ICHX.eq.2) THEN
992	ws=vitesse2(i,k,0)
993	if(i.lt.nz) wi=vitesse2(i+1,k,0)
994	if(i.eq.nz) wi=vitesse2(i ,k,1)
995	w=(ws+wi)/2.
996	zni =zb(i)-w*dt
997	znip1=zb(i)-dzb(i)-w*dt
998	ENDIF
999
1000	explim=30.
1001
1002	IF(ICHX.eq.1 .or. ICHX.eq.2) THEN
1003	ws=vitesse2(i,k,0)
1004	zs=zb(i)
1005	wi=vitesse2(i,k,1)
1006	zi=z(i)
1007
1008	c if(wi.eq.ws) wi=ws/1.001 ! sinon ca plante !
1009	if(abs(wi-ws)/wi .le. 1.e-3) wi=ws/1.001 ! sinon ca plante !
1010
1011	if(wi.ne.0.) alpha= alog(ws/wi) /(zs-zi)
1012	argexp=alpha*zs
1013	if(argexp.lt.-explim) argexp=-explim
1014	if(argexp.gt. explim) argexp=+explim
1015	v0 = ws/exp(argexp)
1016
1017	argexp=alpha*zb(i)
1018	if(argexp.lt.-explim) argexp=-explim
1019	if(argexp.gt. explim) argexp=+explim
1020	arg1=1.+v0alphaexp(argexp)*dt
1021
1022	argexp=alpha*(zb(i)-dzb(i))
1023	if(argexp.lt.-explim) argexp=-explim
1024	if(argexp.gt. explim) argexp=+explim
1025	arg2=1.+v0alphaexp(argexp)*dt
1026
1027	iiter=0
1028	101 continue
1029
1030	if (iiter.le.25) then
1031	if(arg1.le.0. .or. arg2.le.0.) then
1032	print*,ihor,i,iiter, 'ajustement vitesse',arg1,arg2
1033	print*,ws,wi, ' w1 w2 anc valeurs'
1034	print*,alpha, ' alpha anc valeurs'
1035	print*,rad(i,k), 'r(j,k)'
1036	print*,mmu(i,k), ' mmu(j,k)'
1037	print*,t(i), ' t(j,k)'
1038	print*,arg1,arg2, ' arg1 arg2 anc valeurs'
1039	arg2=1+(arg2-1.)/2.
1040	arg1=1+(arg1-1.)/2.
1041	iiter=iiter+1
1042	print*,arg1,arg2, ' arg1 arg2 nle valeurs'
1043	goto 101
1044	endif
1045	else
1046	stop
1047	endif
1048
1049	if(arg1.lt.0.) print*,'arg1:',arg1
1050	if(arg2.lt.0.) print*,'arg2:',arg2
1051
1052	deltazs=-alog(arg1)/alpha
1053	deltazi=-alog(arg2)/alpha
1054
1055	zni =zb(i)+deltazs
1056	znip1=zb(i)-dzb(i)+deltazi
1057
1058	ENDIF
1059
1060	if(zni.ne.znip1) xft=zni/(zni-znip1)
1061	if(zni.eq.znip1 .and. zni.le.0.) xft=0.
1062	if(zni.eq.znip1 .and. zni.gt.0.) then
1063	xft=0.
1064	print*,'zni..znip1', zni,znip1
1065	endif
1066
1067	* Si des aerosols touchent le sol (zni < 0 ) alors on fixe
1068	* le niveau a 0, et on elimine les aerosols correspondants
1069	*-----------------------------------------------------------
1070
1071	if(znip1 .lt. 0.) znip1=0.
1072	if(zni .lt. 0.) zni =0.
1073
1074	if(zni.le.0. .and. znip1 .le. 0.) then
1075	c print*,'voie 1 / disparition complete'
1076	xft=0.
1077	xf=1.
1078	xcnt=xcnt+xf
1079	puit(i)=puit(i)+xf
1080	endif
1081
1082	if(zni.gt.0. .and. znip1 .le. 0.) then
1083	c print*,'voie 2 / disparitipon partielle'
1084	xf=(1-xft)
1085	xcnt=xcnt+xf
1086	puit(i)=puit(i)+xf
1087	endif
1088
1089	if(zni.gt.0. .and. znip1 .gt. 0.) then
1090	c print*,'voie 3 / pas de disparition'
1091	xft=1.
1092	xf=(1-xft)
1093	xcnt=xcnt+xf
1094	puit(i)=puit(i)+xf
1095	endif
1096
1097	jsup=nz+1
1098	jinf=nz+1
1099	do j=nztop,nz
1100	if(zni.le.zb(j) .and. zni.ge.zb(j)-dzb(j)) jsup=j
1101	if(znip1.le.zb(j).and. znip1.ge.zb(j)-dzb(j)) jinf=j
1102	enddo
1103	if(zni .ge. 0. .and. zni .lt. 1.e-3) jsup=nz
1104	if(znip1 .ge. 0. .and. znip1 .lt. 1.e-3) jinf=nz
1105
1106
1107	* Volume inclu dans un seul niveau
1108	*----------------------------------
1109
1110	if (jsup .eq. jinf. and. jsup.ge. nz+1) then
1111	xcnt=xcnt+1.
1112	print*,'cas impossible'
1113	print*,'alpha= ',alpha
1114	print*,'ws wi ',ws,wi
1115	print*,'deltazs deltazi ',deltazs,deltazi
1116	print*,' r(i,k) mmu(i,k)', rad(i,k),mmu(i,k)
1117	print*,' t(j,k)',t(i)
1118	print*,zni,znip1,jsup,jinf
1119	print*,'STOP'
1120	STOP
1121	endif
1122
1123	if (jsup .eq. jinf. and. jsup.le. nz) then
1124	xf=1.
1125	xcnt=xcnt+xft*xf
1126	if(jinf.le.nz) then
1127	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1128	endif
1129	endif
1130
1131	* Volume a cheval sur 2 niveaux
1132	*--------------------------------
1133
1134	if (jinf .eq. jsup+1) then
1135	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
1136	xcnt=xcnt+xf*xft
1137	if(jsup.le.nz) then
1138	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
1139	endif
1140	xf=(zb(jinf)-znip1)/(zni-znip1)
1141	xcnt=xcnt+xf*xft
1142	if(jinf.le.nz) then
1143	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1144	endif
1145	endif
1146
1147	* Volume a cheval sur 3 ou plus de niveaux
1148	*------------------------------------------
1149
1150	if (jinf .gt. jsup+1) then
1151
1152	c print*,' voie C / dans N cases'
1153	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
1154	xcnt=xcnt+xf*xft
1155	if(jsup.le.nz) then
1156	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
1157	endif
1158
1159	xf=(zb(jinf)-znip1)/(zni-znip1)
1160	xcnt=xcnt+xf*xft
1161	if(jinf.le.nz) then
1162	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
1163	endif
1164
1165	do jj=jsup+1,jinf-1
1166	xf=(zb(jj)-zb(jj+1))/(zni-znip1)
1167	xcnt=xcnt+xf*xft
1168	if(jj.le.nz) then
1169	echange(jj,i,ihor)=echange(jj,i,ihor)+xft*xf
1170	endif
1171	enddo
1172
1173	endif
1174
1175
1176	* et sur les rayons...
1177	*---------------------
1178	30 continue
1179
1180	* fin de la grande boucle sur les altitudes...
1181	*----------------------------------------------
1182	20 continue
1183
1184	* Calcul etat final Cfinal = Echange*initial
1185	*----------------------------------------------
1186
1187	compte=0.
1188	compte2=0.
1189	do j=1,nz
1190	xepl=0.
1191	do jj=1,nz
1192	xepl=xepl+echange(jj,j,ihor)
1193	compte=compte+echange(jj,j,ihor)
1194	compte2=compte2+echange(jj,j,ihor)
1195	enddo
1196	compte2=compte2+puit(j)
1197	enddo
1198
1199	if(abs(compte2-nz) .gt. 1.e-4)
1200	& print*,'Matrice calculee#',ihor,'tx d expl (/55):',
1201	& compte,compte2
1202
1203
1204	* Fin du calcul de la matrice d*echange
1205	*----------------------------------------------
1206
1207	do j=nztop,nz
1208	do k=1,nrad
1209
1210	do jj=nztop,nz
1211	ci(j,k,lf)=ci(j,k,lf)+
1212	& echange(j,jj,ihor)*ci(jj,k,li)
1213	ci1(j,k,lf)=ci1(j,k,lf)+
1214	& echange(j,jj,ihor)*ci1(jj,k,li)
1215	ci2(j,k,lf)=ci2(j,k,lf)+
1216	& echange(j,jj,ihor)*ci2(jj,k,li)
1217	ci3(j,k,lf)=ci3(j,k,lf)+
1218	& echange(j,jj,ihor)*ci3(jj,k,li)
1219	enddo
1220	enddo
1221	enddo
1222
1223
1224	* Controles et affichage Bilan
1225	*----------------------------------------------
1226
1227	bilan11=0.
1228	bilan12=0.
1229	bilan13=0.
1230	bilan14=0.
1231	bilan15=0.
1232	do k=1,nrad
1233	do j=nztop,nz
1234	c(j,k,lf) =ci(j,k,lf)/dzb(j)
1235	c1(j,k,lf)=ci1(j,k,lf)/dzb(j)
1236	c2(j,k,lf)=ci2(j,k,lf)/dzb(j)
1237	c3(j,k,lf)=ci3(j,k,lf)/dzb(j)
1238	bilan15=bilan15+ ci(j,k,lf)
1239	bilan11=bilan11+ci1(j,k,lf)
1240	bilan12=bilan12+ci2(j,k,lf)
1241	bilan13=bilan13+ci3(j,k,lf)
1242	bilan14=bilan14+
1243	& ci(j,k,lf)4./3.pirf(k)3.vrat_e**(k-imono)
1244	enddo
1245	enddo
1246
1247	c print*,'sedifn_fast Bilans:'
1248	c & ,bilan11,bilan12,bilan13
1249	c print*,'Bilan1:',bilan1,bilan11
1250	c print*,'Bilan2:',bilan2,bilan12
1251	c print*,'Bilan3:',bilan3,bilan13
1252	c print*,'Bilan5:',bilan5,bilan15
1253
1254	dice1=0.
1255	dice2=0.
1256	dice3=0.
1257	dice4=0.
1258	dice1=(bilan11-bilan1)rhoi_ch4 !glace 1 m^3.m^-2 kg.m^-3 pourquoi ????
1259	dice2=(bilan12-bilan2)*rhoi_c2h6 !glace 2
1260	dice3=(bilan13-bilan3)*rhoi_c2h2 !glace 3
1261	dice4=(bilan14-bilan4)*rhol !noyaux
1262
1263	return
1264	end
1265
1266
1267	subroutine coagul2(ihor)
1268
1269	*********************************************************
1270	* ce programme calcule la nouvelle concentration dans *
1271	* le a ieme intervalle de rayon, a l'altitude h, a *
1272	* l'instant t+dt *
1273	*********************************************************
1274
1275	use dimphy
1276	IMPLICIT NONE
1277	#include "dimensions.h"
1278	#include "microtab.h"
1279	#include "varmuphy.h"
1280
1281	* declaration des blocs communs
1282	*------------------------------
1283
1284	common/x2ctps/li,lf,dt
1285	common/x2con/c,c1,c2,c3,caer
1286	common/cldpart/rad,mmu
1287
1288	* declaration des variables
1289	* --------------------------
1290
1291	integer li,lf
1292	real dt
1293	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1294	& c3(nz,nrad,2)
1295	real caer(nz,nrad,2)
1296	real rad(nz,nrad),mmu(nz,nrad)
1297
1298	* declaration des variables propres au ss-programme
1299	* -------------------------------------------------
1300
1301	integer h,a,ihor,i
1302	real pr,pe,eta1,eta2
1303	real sum11,sum12,sum13,sum21,sum22,sum23
1304	real rfb,rfx,rpr
1305
1306	* traitement
1307	* ----------
1308
1309
1310	sum11=0.
1311	sum12=0.
1312	sum13=0.
1313	sum21=0.
1314	sum22=0.
1315	sum23=0.
1316
1317	do 10 h=nztop,nz
1318	do 11 a=1,nrad ! boucle aerosol secs
1319	call pertpro2(ihor,h,a,pe,pr)
1320
1321	if(1+dtpe .gt. 1.1) print,a,1+dt*pe,' scav'
1322	caer(h,a,lf)=(caer(h,a,li)+prdt)/(1+dtpe)
1323	c(h,a,lf)=c(h,a,li)
1324	c1(h,a,lf)=c1(h,a,li)
1325	c2(h,a,lf)=c2(h,a,li)
1326	c
1327	c eta1=0.
1328	c eta2=0.
1329	c if(c(h,a,li) .ne. 0.) eta1=c1(h,a,li)/c(h,a,li)
1330	c if(c(h,a,li) .ne. 0.) eta2=c2(h,a,li)/c(h,a,li)
1331	c c(h,a,lf) =( c(h,a,li)+prdt)/(1+dtpe)
1332	c c1(h,a,lf)=(c1(h,a,li)+eta1prdt)/(1+dt*pe)
1333	c c2(h,a,lf)=(c2(h,a,li)+eta2prdt)/(1+dt*pe)
1334	c
1335	c
1336	c sum11=sum11+c(h,a,li)
1337	c sum12=sum12+c1(h,a,li)
1338	c sum13=sum13+c2(h,a,li)
1339	c
1340	c sum21=sum21+ c(h,a,lf)
1341	c sum22=sum22+c1(h,a,lf)
1342	c sum23=sum23+c2(h,a,lf)
1343	c
1344	c
1345	11 continue
1346	10 continue
1347
1348
1349	if (nztop.ne.1) then
1350	do 12 h=1,nztop-1
1351	do 12 a=1,nrad
1352	c(h,a,lf)=c(h,a,li)
1353	c1(h,a,lf)=c1(h,a,li)
1354	c2(h,a,lf)=c2(h,a,li)
1355	caer(h,a,lf)=caer(h,a,li)
1356	12 continue
1357	endif
1358
1359	return
1360	end
1361
1362
1363	*__________________________________________________________________________
1364
1365	subroutine calcoag2
1366
1367	***************************************************************
1368	* *
1369	* Ce programme calcule les coefficients de collection d'une *
1370	* particule de rayon x avec une particule de rayon b a une *
1371	* altitude donnee h *
1372	***************************************************************
1373
1374	* declaration des blocs communs
1375	*------------------------------
1376	use dimphy
1377	IMPLICIT NONE
1378	#include "dimensions.h"
1379	#include "microtab.h"
1380	#include "varmuphy.h"
1381
1382	common/x2ctps/li,lf,dt
1383	common/x2con/c,c1,c2,c3,caer
1384	common/cldpart/rad,mmu
1385	common/x2coag/k
1386	common/x2frac/rfg,dfg
1387
1388	* declaration des variables
1389	* --------------------------
1390
1391	integer li,lf
1392	real dt
1393	real knu2,nud2,k(nz,nrad,nrad)
1394	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1395	& c3(nz,nrad,2)
1396	real caer(nz,nrad,2)
1397	real rad(nz,nrad),mmu(nz,nrad)
1398	real rfg(nz),dfg(nz,nrad)
1399
1400
1401	* declaration des variables propres au ss-programme
1402	* -------------------------------------------------
1403
1404	integer h,b,x,ihor,i
1405	real nua,lambb,lambx,knb,knx,alphab,alphax,d,e,f,kcg
1406	real db,dx,rm,dm,deltab,deltax,del,g,beta,gx,gb
1407	real rfb,rfx,rpr
1408	real*8 ne,qe,epso
1409	real*8 corelec,yy
1410
1411	real kco,vx,vb,vitesse2,sto,ee,a,dd,bb,p0,t0,l0,ccol
1412	real st(37),ef(37)
1413	real vitesse,knu
1414	external vitesse,knu
1415
1416
1417	* initialisation
1418	* --------------
1419
1420
1421
1422	* -nombres de STOCKES
1423
1424	data(st(i),i=1,37)/1.35,1.5,1.65,1.85,2.05,2.25,2.5,2.8,3.1,
1425	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,
1426	s 10.9,11.1,13.5,15.3,17.25,20.5,24.5,30.4,39.3,48,57,86.,
1427	s 187.,600./
1428
1429	* -coef. d'efficacite de collection
1430
1431	ef(1)=3.75
1432	ef(2)=8.75
1433	do 11 i=3,37
1434	ef(i)=ef(i-1)+2.5
1435	11 continue
1436
1437	do 2 i=1,37
1438	ef(i)=ef(i)*1e-2
1439	2 continue
1440
1441	qe=1.6e-19
1442	ne=-30e+6
1443	epso=1e-9/(36*pi)
1444
1445	d=1.257
1446	e=0.4
1447	f=-1.1
1448
1449	* iteration sur z
1450
1451	do 1 h=1,nz
1452
1453	nua=nud2(h,1)
1454
1455	* iteration sur les rayons
1456
1457	do 1 b=1,nrad ! boucle aerosols secs : indice ''
1458
1459	knb=knu(h,b,1) ! knu et vitesse se trouvent dans brume.F
1460	vb=vitesse(h,b,1) ! ils concernent les aerosols
1461
1462	do 1 x=1,nrad !boucles gouttes : indice '2'
1463
1464	knx=knu2(h,x,1)
1465	vx=vitesse2(h,x,1)
1466
1467	COAGULATION **************************************************
1468	--------------**************************************************
1469	* calcul du terme correcteur 'slip-flow'
1470
1471	alphab=d+e*exp(f/knb)
1472	alphax=d+e*exp(f/knx)
1473
1474	* calcul du coefficient de diffusion
1475
1476
1477	rfb=(r_e(b)*(3./df(b))) ((rf(b)) **(1.-3./df(b)))
1478	rfx=(rad(h,x)**
1479	& (3./dfg(h,x)))((rfg(x))*(1.-3./dfg(h,x)))
1480
1481	db=kbzt(h)(1+alphabknb)/(6pinuarfb)
1482	dx=kbzt(h)(1+alphaxknx)/(6pinuarfx)
1483
1484	* calcul du coefficient de coagulation
1485
1486	rpr=rfb+rfx
1487	kcg=4pirpr*(db+dx)
1488
1489	* calcul de la vitesse thermique
1490
1491	gb=sqrt(6kbzt(h)/(rholpi2r_e(b)**3))
1492	gx=sqrt(6kbzt(h)/(rholpi2rad(h,x)**3))
1493
1494	* calcul du libre parcours apparent des aerosols
1495
1496	lambb=8db/(pigb)
1497	lambx=8dx/(pigx)
1498
1499	*calcul du terme correcteur beta
1500
1501	rm=rpr/2.
1502	dm=(dx+db)/2.
1503	g=sqrt(gx2+gb2)
1504	deltab=(((2rfb+lambb)3-(4rfb2+lambb2)**1.5)
1505	s /(6rfblambb)-2rfb)sqrt(2.)
1506	deltax=(((2rfx+lambx)3-(4rfx2+lambx2)**1.5)
1507	s /(6rfxlambx)-2rfx)sqrt(2.)
1508	del=sqrt(deltab2+deltax2)
1509	beta=1/((rm/(rm+del/2))+(4dm/(grm)))
1510
1511	* calcul du coefficient de coagulation corrige
1512
1513	kcg=kcg*beta
1514
1515	COALESCENCE ************************************************
1516	-------------************************************************
1517
1518	kco=0.
1519
1520	if ( b.eq. x) continue ! goto 9
1521
1522
1523	* calcul du nombre de Stockes de la petite particule
1524
1525	sto=2rholrfx*2abs(vx-vb)/(9nuarfb)
1526
1527	* calcul du coef. de Cunningham-Millikan
1528
1529	a=1.246
1530	bb=0.42
1531	dd=0.87
1532	l0=0.653e-7
1533	p0=101325.
1534	t0=288.
1535
1536	ee=1+
1537	& (l0t(h)p0(a+bbexp(-ddrfxt0p(h)/(l0t(h)*p0))))
1538	s /(rfxt0p(h))
1539
1540	* calcul du nombre de Stockes corrige
1541
1542	sto=sto*ee
1543
1544	if (sto .le. 1.2) goto 9
1545
1546	if (sto .ge. 600.) then
1547	ccol=1.
1548	goto 8
1549	endif
1550
1551	* recherche du coefficient de collection
1552
1553	do 3 i=1,37
1554	if (sto .gt. st(i)) then
1555	goto 3
1556	endif
1557	if (sto .eq. st(i)) then
1558	ccol=ef(i+1)
1559	else
1560	ccol=ef(i)
1561	endif
1562	goto 8
1563	3 continue
1564
1565	* calcul du coefficient de coalescence
1566
1567	8 kco=pi(rfb+rfx)2ccol*abs(vb-vx)
1568
1569	9 continue
1570
1571	CORRECTION ELECTRICITE *****************************
1572	------------------------****************************
1573
1574	c yy=1.d0ne2r(x)r(b)qe**2
1575	c & /(1.d0kbzt(h)(r(b)+r(x))4piepso)
1576	c corelec=0.
1577	c if (yy.lt.50.) corelec=yy/(exp(yy)-1.)
1578
1579	corelec=1.
1580
1581	c b=aerosol
1582	c x=gouttes
1583
1584	k(h,b,x)=(kcg+kco)*corelec
1585
1586	c
1587	c ATTENTION, IL N'Y A PLUS DE SYMETRIE...
1588	c k(ihor,h,x,b)=k(ihor,h,b,x)
1589	c
1590	c
1591
1592
1593	1 continue
1594	return
1595	end
1596
1597	*______________________________________________________________________
1598
1599	subroutine pertpro2(ihor,h,a,l_,pr_)
1600
1601	*****************************************************************************
1602	* *
1603	* ce programme permet le calcul du terme de production (pr) et de perte (l)*
1604	* pour le phenomene de coagulation *
1605	* dans le a ieme intervalle de rayon a une altitude h *
1606	****************************************************************************
1607
1608	* declaration des blocs communs
1609	*------------------------------
1610	use dimphy
1611	IMPLICIT NONE
1612	#include "dimensions.h"
1613	#include "microtab.h"
1614	#include "varmuphy.h"
1615
1616	common/x2ctps/li,lf,dt
1617	common/x2con/c,c1,c2,c3,caer
1618	common/cldpart/rad,mmu
1619	common/x2coag/k
1620
1621
1622	* declaration des variables
1623	* --------------------------
1624
1625	integer li,lf
1626	real dt
1627	real dr(nrad),dv(nrad)
1628	real k(nz,nrad,nrad)
1629	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
1630	& c3(nz,nrad,2)
1631	real caer(nz,nrad,2)
1632	real rad(nz,nrad),mmu(nz,nrad)
1633
1634
1635	* declaration des variables propres au ss-programme
1636	* -------------------------------------------------
1637
1638	integer h,b,a,x,ihor,i
1639	real*8 pr,ss,s,l
1640	real pr_,l_,vol,del
1641
1642	* traitement
1643	* -----------
1644
1645
1646
1647	* production
1648	*+++++++++++++
1649	s=0.d0
1650	ss=0.d0
1651	pr=0.d0
1652
1653	c Pas de production d'aerosols par scavenging !!!
1654
1655	pr=0.d0
1656
1657	* perte
1658	*- - - - -
1659
1660	l=0.d0
1661
1662
1663	do 10 x=1,nrad ! boucle sur les gouttes
1664	l=l+k(h,a,x)*c(h,x,li)
1665	if(l.ne.0.d0) then
1666	print*,a,x,k(h,a,x),c(h,x,li),l,' : detail coal'
1667	endif
1668	10 continue
1669
1670
1671	#ifdef CRAY
1672	l_=l
1673	pr_=pr
1674	#else
1675	l_=sngl(l)
1676	pr_=sngl(pr)
1677	#endif
1678
1679	return
1680
1681	end
1682

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

Download in other formats:

Original Format