Context Navigation

brume3D.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: 30.4 KB

Line
1	subroutine brume(ngrid,tab1,x1,
2	& xnz,xnrad,taused,ihor,
3	& pmu0,pfract,
4	& precip)
5
6
7	*1 c nombre de particules de la grille de rayon r a l'altitude z
8	*2 dt, pas de temps, en heure.
9
10	--------------------------------------------------------------
11	* *
12	* ENTRE 0 ET 1000 KILOMETRES *
13	* *
14	* la dimension fractale est en tableau, attention au *
15	* raccordement entre le regime moleculaire et le regime *
16	* fluide *
17	* *
18	* Modele microphysique: Cabane et al.,1992 / *
19	* Modele version fractale: Cabane et al.,1993 / *
20	* *
21	--------------------------------------------------------------
22	* VERSION DU 2 JUIN 1993 --- AUT 1994 --- 11/04/96
23	*
24	* changer: altitude de production z0=/taux de production ctot=
25	* : la charge/micron, ne
26	* : df(h),rf...
27	* raccordement aknc
28	*
29	* declaration des blocs communs
30	*------------------------------
31
32	use dimphy
33	IMPLICIT NONE
34	#include "dimensions.h"
35	#include "microtab.h"
36	#include "varmuphy.h"
37
38
39	real pmu0,pfract
40
41	common/ctps/li,lf,dt
42	common/con/c
43	common/coag/k
44	common/effets/xsaison
45
46
47	* declaration des variables communes
48	* ----------------------------------
49
50	integer xnz,xnrad,ngrid
51	integer li,lf,ihor
52	real dt
53	real c0(nz,nrad),c(nz,nrad,2)
54	real k(nz,nrad,nrad),knu
55	real xsaison
56	real taused(nz,nrad)
57	real precip(ngrid,5)
58
59
60
61	* variables internes
62	* ------------------
63
64	integer h,ti,itime,i,j
65	real tab1(nz,nrad)
66	real x1
67	real somme,v1,dice3
68
69	real vitesse
70
71	save itime
72	data itime/0/
73
74	* initialisation
75	* --------------
76
77	* effet saisonnier
78	* ----------------
79
80
81	xsaison=0.
82	xsaison=pmu04.pfract
83	!=Pi si fract=1/2 (equinoxe) et
84	! si mu0(ihor)=1 sous le soleil
85	! exactement.
86
87	c xsaison=0.
88	c if (ihor.le.9.or.ihor.ge.41) xsaison=8. ! rapport des surfaces
89	c xsaison=1.
90
91	* controles
92	* ---------
93
94	if (nrad.ne.xnrad) stop 'nrad.ne.xnrad'
95	if (nz.ne.xnz) stop 'nz.ne.xnz'
96
97	do i=1,nz
98	do j=1,nrad
99	c(i,j,1)=tab1(i,j)
100	c(i,j,2)=0.0
101	enddo
102	enddo
103
104	dt=x1
105
106	if (itime.eq.0) then
107	ITIME=1
108	call init
109	call calcoag
110	endif
111
112	do i=1,nz,1
113	do j=1,nrad
114	v1=vitesse(i,j,0)
115	taused(i,j)=(8.314*t(i)/28.e-3/1.35)/v1
116	enddo
117	enddo
118
119	call coagul
120
121	call production(ihor)
122
123	li=3-li
124	lf=3-lf
125
126	call sedif(dice3)
127	c En theorie, dice3 est NEGATIF (en sedimentant on ne fait que perdre des aerosols)
128	c Le precipitations sont comptées positivement.
129	precip(ihor,5)=-dice3/rhol ! m3/m2=m
130
131	li=3-li
132	lf=3-lf
133
134	do i=1,nz
135	do j=1,nrad
136	tab1(i,j)=c(i,j,li) ! li=1
137	enddo
138	enddo
139
140	return
141
142	end
143	*________________________________________________________________________
144
145	subroutine coagul
146
147	*********************************************************
148	* ce programme calcule la nouvelle concentration dans *
149	* le a ieme intervalle de rayon, a l'altitude h, a *
150	* l'instant t+dt *
151	*********************************************************
152	use dimphy
153	IMPLICIT NONE
154	#include "dimensions.h"
155	#include "microtab.h"
156	#include "varmuphy.h"
157
158
159	* declaration des blocs communs
160	*------------------------------
161
162	common/ctps/li,lf,dt
163	common/con/c
164
165	* declaration des variables
166	* --------------------------
167
168	integer li,lf
169	real dt
170	real c(nz,nrad,2)
171
172	* declaration des variables propres au ss-programme
173	* -------------------------------------------------
174
175	integer h,a
176	real pr,pe
177
178	* traitement
179	* ----------
180
181	do h=nztop,nz
182	do a=1,nrad
183	call pertpro(h,a,pe,pr)
184	c if((1+dt*pe).lt.0.) stop 'denom.eq.0'
185	c(h,a,lf)=(c(h,a,li)+prdt)/(1+dtpe)
186	enddo
187	enddo
188
189	if (nztop.ne.1) then
190	do h=1,nztop-1
191	do a=1,nrad
192	c(h,a,lf)=c(h,a,li)
193	enddo
194	enddo
195	endif
196
197	return
198	end
199
200
201	*__________________________________________________________________________
202
203	subroutine calcoag
204
205	***************************************************************
206	* *
207	* Ce programme calcule les coefficients de collection d'une *
208	* particule de rayon x avec une particule de rayon b a une *
209	* altitude donnee h *
210	***************************************************************
211
212	* declaration des blocs communs
213	*------------------------------
214	use dimphy
215	IMPLICIT NONE
216	#include "dimensions.h"
217	#include "microtab.h"
218	#include "varmuphy.h"
219
220	common/ctps/li,lf,dt
221	common/con/c
222	common/coag/k
223
224	* declaration des variables
225	* --------------------------
226
227	integer li,lf,i
228	real dt
229	real c(nz,nrad,2)
230	real knu,nud,k(nz,nrad,nrad)
231
232	* declaration des variables propres au ss-programme
233	* -------------------------------------------------
234
235	integer h,b,x
236	real nua,lambb,lambx,knb,knx,alphab,alphax,d,e,f,kcg
237	real db,dx,rm,dm,deltab,deltax,del,g,beta,gx,gb
238	real rfx,rfb,rpr
239	real*8 ne,qe,epso
240	real*8 corelec,yy
241
242	real kco,vx,vb,vitesse,sto,ee,a,dd,bb,p0,t0,l0,ccol
243	real st(37),ef(37)
244
245	* initialisation
246	* --------------
247	c print,'*** calcoag'
248
249	* -nombres de STOCKES
250
251	data(st(i),i=1,37)/1.35,1.5,1.65,1.85,2.05,2.25,2.5,2.8,3.1,
252	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,
253	s 10.9,11.1,13.5,15.3,17.25,20.5,24.5,30.4,39.3,48,57,86.,
254	s 187.,600./
255
256	* -coef. d'efficacite de collection
257
258	ef(1)=3.75
259	ef(2)=8.75
260	do i=3,37
261	ef(i)=ef(i-1)+2.5
262	enddo
263
264	do i=1,37
265	ef(i)=ef(i)*1e-2
266	enddo
267
268	qe=1.6e-19
269	ne=-30.e+6 ! "vieille valeur!"
270	ne=-15.e+6 ! pour fitter DISR !
271
272	epso=1e-9/(36*pi)
273
274	d=1.257
275	e=0.4
276	f=-1.1
277
278	* iteration sur z
279
280	do 1 h=1,nz
281	nua=nud(h,1)
282
283	* iteration sur les rayons
284
285	do 1 b=1,nrad
286
287	knb=knu(h,b,1)
288	vb=vitesse(h,b,1)
289
290	do 1 x=1,b
291
292	knx=knu(h,x,1)
293	vx=vitesse(h,x,1)
294
295	COAGULATION **************************************************
296	--------------**************************************************
297	* calcul du terme correcteur 'slip-flow'
298
299	alphab=d+e*exp(f/knb)
300	alphax=d+e*exp(f/knx)
301
302	* calcul du coefficient de diffusion
303
304	rfb=(r_e(b)*(3./df(b)))((rf(b))**(1.-3./df(b)))
305	rfx=(r_e(x)*(3./df(x)))((rf(x))**(1.-3./df(x)))
306	db=kbzt(h)(1+alphabknb)/(6pinuarfb)
307	dx=kbzt(h)(1+alphaxknx)/(6pinuarfx)
308
309	* calcul du coefficient de coagulation
310
311	rpr=rfb+rfx
312	kcg=4pirpr*(db+dx)
313
314	* calcul de la vitesse thermique
315
316	gx=sqrt(6kbzt(h)/(rholpi2r_e(x)**3))
317	gb=sqrt(6kbzt(h)/(rholpi2r_e(b)**3))
318
319	* calcul du libre parcours apparent des aerosols
320
321	lambb=8db/(pigb)
322	lambx=8dx/(pigx)
323
324	*calcul du terme correcteur beta
325
326	rm=rpr/2.
327	dm=(dx+db)/2.
328	g=sqrt(gx2+gb2)
329	deltab=(((2rfb+lambb)3-(4rfb2+lambb2)**1.5)
330	s /(6rfblambb)-2rfb)sqrt(2.)
331	deltax=(((2rfx+lambx)3-(4rfx2+lambx2)**1.5)
332	s /(6rfxlambx)-2rfx)sqrt(2.)
333	del=sqrt(deltab2+deltax2)
334	beta=1/((rm/(rm+del/2))+(4dm/(grm)))
335
336	* calcul du coefficient de coagulation corrige
337
338	kcg=kcg*beta
339	c print*,' kcg:', rfb,rfx,knb,knx,db,dx
340	c print*,' beta:', gx,gb,lambb,lambx,rm,dm
341	c print*,' beta:', g,deltab,deltax,del,beta
342	c print*,' '
343
344	COALESCENCE ************************************************
345	-------------************************************************
346
347	kco=0.
348
349	if (b.eq.x) goto 9
350
351	* calcul du nombre de Stockes de la petite particule
352
353	sto=2rholrfx*2abs(vx-vb)/(9nuarfb)
354
355	* calcul du coef. de Cunningham-Millikan
356
357	a=1.246
358	bb=0.42
359	dd=0.87
360	l0=0.653e-7
361	p0=101325.
362	t0=288.
363
364	ee=1+(l0t(h)p0*
365	& (a+bbexp(-ddrfxt0p(h)/(l0t(h)p0))))
366	s /(rfxt0p(h))
367
368	* calcul du nombre de Stockes corrige
369
370	sto=sto*ee
371
372	if (sto .le. 1.2) goto 9
373
374	if (sto .ge. 600.) then
375	ccol=1.
376	goto 8
377	endif
378
379	* recherche du coefficient de collection
380
381	do 3 i=1,37
382	if (sto .gt. st(i)) then
383	goto 3
384	endif
385	if (sto .eq. st(i)) then
386	ccol=ef(i+1)
387	else
388	ccol=ef(i)
389	endif
390	goto 8
391	3 continue
392
393	* calcul du coefficient de coalescence
394
395	8 kco=pi(rfb+rfx)2ccol*abs(vb-vx)
396
397	9 continue
398
399	CORRECTION ELECTRICITE *****************************
400	------------------------****************************
401
402	yy=1.d0ne2r_e(x)r_e(b)qe**2
403	& /(1.d0kbzt(h)(r_e(b)+r_e(x))4piepso)
404
405
406	corelec=0.
407	if (yy.lt.50.) corelec=yy/(exp(yy)-1.)
408	if (yy.le.1.e-3) corelec=1.
409
410	k(h,b,x)=(kcg+kco)*corelec
411	k(h,x,b)=k(h,b,x)
412
413
414	1 continue
415	return
416	end
417
418	*______________________________________________________________________
419
420	real function lambda(j,indic)
421	*
422	------------------------------------------------------------------
423	* fonction calculant le libre parcours moyen des molecules *
424	* atmospheriques( rayon =ra) se trouvant dans la couche no j. *
425	* pour indic=0 ...... la particule se trouve a la frontiere entre*
426	* les couches j et j-1 *
427	* pour indic=1 ...... la particule se trouve au milieu de la *
428	* la couche j *
429	------------------------------------------------------------------
430	*
431	* declaration des blocs communs
432	*------------------------------
433	IMPLICIT NONE
434	#include "dimensions.h"
435	#include "microtab.h"
436	#include "varmuphy.h"
437
438	* declaration des variables communes
439	* ----------------------------------
440
441	integer i,j
442
443	* declaration des variables internes
444	* ----------------------------------
445
446	integer indic
447	real pp,ra
448
449	ra=1.75e-10
450
451	* traitement
452	* ----------
453
454	if (indic.eq.0) then
455	pp=pb(j)
456	else
457	if (indic.ne.1) then
458	print*,'erreur argument fonction lambda'
459	return
460	endif
461	pp=p(j)
462	endif
463
464	lambda=kbzt(j)/(4sqrt(2.)pi(ra*2)pp)
465	end
466
467	*******************************************************************************
468
469	real function knu(j,k,indic)
470	*
471	--------------------------------------------------------------
472	* fonction calculant le nombre de knudsen d'une particule *
473	* d'aerosol de rayon r_e(k) se trouvant dans la couche no j *
474	* indic ...... idem function lambda *
475	--------------------------------------------------------------
476	*
477	* declaration des blocs communs
478	*------------------------------
479	IMPLICIT NONE
480	#include "dimensions.h"
481	#include "microtab.h"
482	#include "varmuphy.h"
483
484	* declaration des variables internes
485	* ----------------------------------
486
487	integer indic,j,k
488	real lambda,rfk
489
490	* traitement
491	* ----------
492
493	if (indic.ne.0 .and.indic.ne.1) then
494	print*,'erreur argument fonction knu'
495	return
496	endif
497
498	rfk=(r_e(k)*(3./df(k)))((rf(k))**(1.-3./df(k)))
499	knu=lambda(j,indic)/rfk
500	end
501
502	*****************************************************************************
503
504	real function nud(j,indic)
505	*
506	--------------------------------------------------------------
507	* fonction calculant la viscosite dynamique (en USI) de l'air *
508	* d'apres la formule de Sutherlant a l'altitude j *
509	* indic ......... idem fonction lambda *
510	--------------------------------------------------------------
511	*
512	IMPLICIT NONE
513	#include "dimensions.h"
514	#include "microtab.h"
515	#include "varmuphy.h"
516
517
518	integer indic,j
519	real nud0,c,tt
520	*
521	nud0=1.74e-5
522	c=109.
523
524	if(indic.ne.0.and.indic.ne.1) then
525	print*,'erreur argument fonction nud'
526	return
527	endif
528
529	if(indic.eq.0) tt=tb(j)
530	if (indic.eq.1) tt=t(j)
531
532	nud=nud0sqrt(tt/293)(1+c/293)/(1+c/tt)
533	end
534
535	****************************************************************************
536
537	real function vitesse(j,k,indic)
538	*
539	-----------------------------------------------------------------
540	* fonction calculant la vitesse de chute d'une particule de rayon*
541	* k se trouvant a l'altitude j suivant la valeur du nombre de *
542	* Knudsen *
543	* indic ....... idem function lambda *
544	-----------------------------------------------------------------
545	*
546	* declaration des blocs communs
547	*------------------------------
548	IMPLICIT NONE
549	#include "dimensions.h"
550	#include "microtab.h"
551	#include "varmuphy.h"
552
553	* declaration des variables internes
554	* ----------------------------------
555
556	integer indic,j,k
557	real w,g,m,a0,zz,nud,knud,tt,rhoh
558	real rbis, rfk,vb,zbx
559	real akncx,knu
560
561	* traitement
562	* ----------
563
564	if (indic.ne.0.and.indic.ne.1) then
565	print*,'erreur argument fonction vitesse'
566	return
567	endif
568
569	if(indic.eq.0) then
570	zz=z(j)+dz(j)/2.
571	tt=tb(j)
572	rhoh=rhob(j)
573	endif
574	if(indic.eq.1) then
575	zz=z(j)
576	tt=t(j)
577	rhoh=rho(j)
578	endif
579
580	g=g0(rtit/(rtit+zz))*2
581	a0=0.74
582	m=(ach4(j)mch4+aar(j)mar+an2(j)*mn2)/nav
583	knud=knu(j,k,indic)
584
585
586	rfk=(r_e(k)*(3./df(k)))((rf(k))**(1.-3./df(k)))
587	c rfk=r_e(7)
588
589
590	w=2./9.rfk(df(k)-1.)rf(k)*(3.-df(k))g*rhol/nud(j,indic)
591
592	w=w(1+1.2517knud+0.4knudexp(-1.1/knud))
593
594
595	c if (p(j).lt.500..and.k.eq.nrad) then
596	c w=0.
597	c endif
598
599	vitesse=w
600
601	end
602	***********************************************************************
603
604	real function kd(h)
605	*
606	--------------------------------------------------------------------
607	* cette fonction calcule le coefficient du terme de 'eddy diffusion'*
608	* a l altitude j *
609	--------------------------------------------------------------------
610	*
611	IMPLICIT NONE
612	#include "dimensions.h"
613	#include "microtab.h"
614	#include "varmuphy.h"
615
616	real zbx
617
618	integer h
619
620	zbx=z(h)+dz(h)/2.
621	if(zbx.le.42000.) then
622	kd=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
623	kd=4.
624	else
625	kd=0.0*kd
626	kd=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
627	endif
628
629	if(zbx.le.50000.) then
630	kd=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
631	endif
632
633	kd=0.0*kd
634
635	return
636	end
637
638
639	*____________________________________________________________________________
640
641	subroutine init
642	*
643	--------------------------------------------------------------------
644	* cette routine effectue : *
645	* 1) interpolation a partir des donnees initiales des *
646	* valeurs de p,t,rho,ach4,aar,an2 sur la grille *
647	* 2) initialisation des constantes (common/phys/) *
648	* 3) initialisation des variables temporelles (common *
649	* /temps/) *
650	* 4) definition des grilles en rayon et verticale *
651	* 5) initialisation de c(z,r,t) avec les donnees du *
652	* fichier unit=1 *
653	* *
654	* les donnees sont des valeurs caracterisques de l'atmosphere de *
655	* TITAN ( voir Lelouch and co ) *
656	--------------------------------------------------------------------
657
658	* declaration des blocs communs
659	*------------------------------
660	use dimphy
661	IMPLICIT NONE
662	#include "dimensions.h"
663	#include "microtab.h"
664	#include "varmuphy.h"
665
666
667	common/ctps/li,lf,dt
668	common/con/c
669
670	* declaration des variables communes
671	* ----------------------------------
672
673	real c(nz,nrad,2)
674	integer li,lf
675	integer i,ii
676	real dt
677
678	* declaration des variables internes
679	* ----------------------------------
680
681	integer nzd
682	parameter (nzd=254)
683	integer limsup,liminf,j1,j2
684	real zd(nzd),ach4d(nzd),rap
685	real m
686
687
688	* initialisation des variables temporelles
689	* ----------------------------------------
690
691	li=1
692	lf=2
693
694
695	* interpolation de xch4,xar et xn2 sur la grille
696	* ----------------------------------------------
697
698	* donnees initiales (Lellouch et al,87)
699	* -------------------------------------
700
701	c print,'***** init'
702	do 1 i=1,168
703	zd(i)=(1000.-5(i-1))1000.
704	1 continue
705	do 2 i=1,78
706	zd(168+i)=(160.-2(i-1))1000.
707	2 continue
708	do 3 i=1,4
709	zd(246+i)=(5.-(i-1))*1000.
710	3 continue
711	do 4 i=1,4
712	zd(250+i)=(1.5-(i-1)0.5)1000.
713	4 continue
714
715	data (ach4d(i),i=1,168)/168*1.5e-2/
716	data (ach4d(i),i=169,254)/63*1.5e-2,1.6e-2,1.8e-2,1.8e-2,
717	& 1.9e-2,2.e-2,2.1e-2,2.3e-2,2.5e-2,2.8e-2,3.1e-2,3.6e-2,
718	& 4.1e-2,4.7e-2,5.7e-2,6.7e-2,7.5e-2,7*8.e-2/
719
720	liminf=0
721	limsup=0
722
723	* interpolation des taux de melange de ch4,ar,n2
724	*-----------------------------------------------
725
726	! do 20 j1=1,nz
727	! do 21 j2=1,nzd
728	! if( zd(j2).le.z(j1)) goto 22
729	!21 continue
730	!22 liminf=j2
731
732	do 20 j1=1,nz
733	do 21 j2=1,nzd
734	if( zd(j2).le.z(j1)) goto 22
735	21 continue
736	22 if (j2.ge.254) j2=254
737	liminf = j2
738
739	if (zd(liminf).eq.z(j1) )then
740	ach4(j1)=ach4d(liminf)
741	goto 20
742	endif
743	if (j2.ne.1) then
744	limsup=j2-1
745	else
746	limsup=j2
747	endif
748
749	if (limsup.eq.liminf) then
750	ach4(j1)=ach4(limsup)
751	else
752	ach4(j1)=ach4d(liminf)-(ach4d(limsup)-ach4d(liminf))/
753	s (zd(limsup)-zd(liminf))*(zd(liminf)-z(j1))
754	endif
755	20 continue
756
757	* rap= aar/an2 cst sur l'altitude
758
759	rap=0.191
760	do 23 i=1,nz
761	an2(i)=(1.-ach4(i))/(1.+rap)
762	aar(i)=rap*an2(i)
763	23 continue
764
765	do 24 i=1,nz
766	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
767	rho(i)=p(i)m/(rgpt(i))
768	24 continue
769
770	do 34 i=1,nz
771	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
772	rhob(i)=pb(i)m/(rgptb(i))
773	c print*,pb(i),m,rgp,tb(i),rhob(i),rho(i)
774	34 continue
775
776	* fin d'interpolation des taux de melange
777	*----------------------------------------
778
779	c print,'*** fin init'
780	540 continue
781	return
782
783	500 print*,'erreur lecture initialisation de c...erreur=',ii
784	stop
785
786	end
787
788	*____________________________________________________________________________
789
790	subroutine pertpro(h,a,l_,pr_)
791
792	*****************************************************************************
793	* *
794	* ce programme permet le calcul du terme de production (pr) et de perte (l)*
795	* pour le phenomene de coagulation *
796	* dans le a ieme intervalle de rayon a une altitude h *
797	*****************************************************************************
798
799	* declaration des blocs communs
800	*------------------------------
801	use dimphy
802	IMPLICIT NONE
803	#include "dimensions.h"
804	#include "microtab.h"
805	#include "varmuphy.h"
806
807
808	common/ctps/li,lf,dt
809	common/con/c
810	common/coag/k
811
812	* declaration des variables
813	* --------------------------
814
815	integer li,lf
816	real dt
817	real c(nz,nrad,2),k(nz,nrad,nrad)
818
819	* declaration des variables propres au ss-programme
820	* -------------------------------------------------
821
822	integer h,b,a,x,i
823	real*8 pr,ss,s,l
824	real pr_,l_,vol,del
825
826	* traitement
827	* -----------
828
829	* production
830	*+++++++++++++
831	s=0.d0
832	ss=0.d0
833	pr=0.
834
835	if (a .eq. 1) goto 2
836	b=a-1
837
838	if (c(h,b,lf) .eq. 0 .and. c(h,b,li) .eq. 0) goto 2
839
840	do 1 i=1,b
841
842	if(c(h,i,li) .eq. 0 .and. c(h,i,lf) .eq. 0) goto 1
843
844	if (i .ne. b)del=1.
845	if (i .eq. b) del=.5
846
847	s=(v_e(i)1.d0)del(k(h,b,i)1.d0)(c(h,i,li)1.d0)+s
848	ss=(v_e(i)1.d0)del(k(h,b,i)1.d0)(c(h,i,lf)1.d0)+ss
849
850	c if (a.eq.2) print*,'SS>',v_e(i),k(h,b,i),c(h,b,lf)
851	c if (a.eq.2) print,'SS>',delv_e(i)k(h,b,i)c(h,b,lf)
852
853
854	1 continue
855
856	* calcul du terme de production
857
858	pr=(c(h,b,lf)s/(vrat_e-1.)+c(h,b,li)ss)/v_e(a)
859	c if (a.eq.2) print*,'PR>',s,ss,c(h,b,lf),v_e(a)
860
861	2 continue
862
863
864	* perte
865	*- - - - -
866
867	l=0
868
869
870	* condition limite : pas de perte dans le dernier intervalle
871
872	if (a .eq. nrad) goto 9
873
874	do 10 x=1,nrad
875
876	if (c(h,x,li) .eq. 0) goto 10
877
878	if (a .lt. x) vol=1.
879	if (a .eq. x) vol=.5*vrat_e/(vrat_e-1)
880	if (a .gt. x) vol=v_e(x)/(v_e(a)*(vrat_e-1))
881
882	l=l+k(h,a,x)c(h,x,li)vol*1.d0
883
884
885	10 continue
886	9 continue
887
888	#ifdef CRAY
889	l_=l
890	pr_=pr
891	#else
892	l_=sngl(l)
893	pr_=sngl(pr)
894	#endif
895	c l_=sngl(l)
896	c pr_=sngl(pr)
897	c if (a.eq.2) print*,'pr_,l_',h,a,pr_,l_
898	c if (a.eq.2) print*,'-----------------------'
899	c if (a.eq.2) STOP
900
901
902
903	return
904
905	end
906
907	*_____________________________________________________________________________
908
909	subroutine production(ihor)
910	*
911	--------------------------------------------------------------------
912	* routine calculant le terme de production des molecules organiques *
913	* composant les aerosols . rini= rayon des aerosols initiaux *
914	--------------------------------------------------------------------
915	*
916	use dimphy
917	IMPLICIT NONE
918	#include "dimensions.h"
919	#include "microtab.h"
920	#include "varmuphy.h"
921	#include "clesphys.h"
922
923	c#include "aerprod.h"
924
925
926	integer ndz,i,ihor,k,k1
927	real c(nz,nrad,2)
928	integer li,lf
929	real dt
930	real zprod,zy,c0,ctot,prod,rini,rfron
931	real xsaison,p0
932	common/ctps/li,lf,dt
933	common /con/c
934	common/effets/xsaison
935
936	! Pressure level of aerosol production
937	p0=p_prodaer
938
939	do i=1,nz-1
940	if (pb(i).lt.p0.and.pb(i+1).gt.p0) zprod=(z(i)+z(i+1))/2.
941	enddo
942
943	ctot=3.5e-13*tx ! ATTENTION, ??COHERENT AVEC INITPAR??
944	ctot=ctot*xsaison !
945
946	c z0=385.e+3
947	zy=20.e+3
948	rini=1.3e-9
949	rini=r_e(1)
950	ndz=50
951	*
952	do 10 i=1,nrad
953	if(rini.lt.r_e(i)) goto 100
954	10 continue
955	100 continue
956	if (i.eq.1) then
957	rini=r_e(1)
958	else
959	rfron=(r_e(i)+r_e(i-1))/2
960	if (rini .lt.rfron) then
961	rini=r_e(i-1)
962	i=i-1
963	else
964	rini=r_e(i)
965	endif
966	endif
967	*
968	c0=ctot/(sqrt(2.pi)zy)
969	c0=c03./(4.pirholrini**3)
970	*
971	do 20 k=nztop,nz
972	prod=0.
973	do 201 k1=1,ndz
974	prod=prod+c0exp(-0.5(((z(k)+dz(k)/2.-k1dz(k)/(2.ndz)
975	s -zprod)/zy)*2))dt/ndz
976	201 continue
977
978	if (prod .le. 1) prod=0.
979	c(k,i,lf)=c(k,i,lf)+prod
980	20 continue
981
982	! do 30 k=nztop,nz
983	! c(k,i,lf)=c(k,i,lf)+prodaer(ihor,nz-k+1,1)
984	! & 3./(4.pirholrini**3)
985	!30 continue
986
987	return
988	end
989
990	-------------------------------------------------------------------
991
992
993
994	subroutine sedif(dice3)
995	*
996	------------------------------------------------------------------
997	* cette routine calcule l'evolution de la fonction de distribution*
998	* c(z,r,t) pour les phenomenes de sedimentation et de diffusion *
999	------------------------------------------------------------------
1000	*
1001	*
1002	* declaration des blocs communs
1003	*------------------------------
1004	use dimphy
1005	IMPLICIT NONE
1006	#include "dimensions.h"
1007	#include "microtab.h"
1008	#include "varmuphy.h"
1009
1010	common/ctps/li,lf,dt
1011	common/con/c
1012
1013	* declaration des variables communes
1014	* ----------------------------------
1015
1016	integer li,lf
1017	integer i,j,k,nb
1018	real dt
1019	real c(nz,nrad,2),dice3,bilan4,bilan14
1020
1021	* declaration des variables internes
1022	* ----------------------------------
1023
1024	real w,w1,dzbX,dc
1025	double precision sigma,theta,hc,l,rap,cmp,wp
1026	double precision fs(nz+1),ft(nz+1)
1027	real as(nz),bs(nz),cs(nz),ds(nz)
1028	double precision asi(nztop:nz),bsi(nztop:nz),
1029	& csi(nztop:nz)
1030	double precision dsi(nztop:nz),xsol(nztop:nz)
1031	real vitesse,kd
1032
1033	external dtridgl
1034	* resolution
1035	*------------
1036
1037
1038	bilan4=0.
1039	do k=1,nrad
1040	do j=nztop,nz
1041	bilan4=bilan4+c(j,k,li)dzb(j)
1042	& 4./3.pirf(k)*3.vrat_e**(k-imono)
1043	enddo
1044	enddo
1045
1046
1047	do 10 k=1,nrad
1048	do 20 j=nztop,nz
1049
1050	if (j.eq.1) goto 20
1051
1052	* calcul de la vitesse corrigee
1053
1054	dzbX=(dz(j)+dz(j-1))/2.
1055	w=-1*vitesse(j,k,0)
1056	if (kd(j).ne.0.) then
1057	theta=0.5(wdzbX/kd(j)+log(rho(j-1)/rho(j)))
1058	if (theta.ne.0) then
1059	sigma=1./dtanh(theta)-1./theta
1060	else
1061	sigma=1.
1062	endif
1063	else
1064	sigma=1.
1065	endif
1066	if(c(j,k,li).eq.0.) then
1067	rap=10.
1068	else
1069	rap=c(j-1,k,li)/c(j,k,li)
1070	if( rap.gt.10.) rap=10.
1071	if( rap.lt.0.1) rap=0.1
1072	endif
1073	if (rap.gt.0.9 .and. rap.lt.1.1) then
1074	w1=w
1075	else
1076	if(w.ne.0) then
1077	hc=dzbX/dlog(rap)
1078	l=dzbX/(wdt)(dexp(-w*dt/hc)-1.)/(1.-rap)
1079	wp=w*1.d0
1080	cmp=dlog(-wp)+abs(sigma)*dlog(l)
1081	if (cmp.gt.38) then
1082	goto 20
1083	endif
1084	w1=-dexp(cmp)
1085
1086	else
1087	w1=0.
1088	endif
1089	endif
1090
1091	c if(w1.ge.0. .or. w1.le.0.) then
1092	c continue
1093	c else
1094	c print*,j,k,w1,hc,dzbx,rap,dlog(rap),' w1'
1095	c endif
1096
1097	* calcul des flux aux interfaces
1098
1099
1100	if (kd(j).ne.0.) then
1101	if (theta.ne.0.) then
1102	ft(j)=(w1+log(rho(j-1)/rho(j))kd(j)/dzbX)/(dexp(2.
1103	s theta)-1.)
1104	fs(j)=ft(j)dexp(2.theta)
1105	else
1106	ft(j)=kd(j)/dzbX
1107	fs(j)=kd(j)/dzbX
1108	endif
1109	else
1110	if (w1.lt.0.)then
1111	ft(j)=-w1
1112	fs(j)=0.
1113	else
1114	ft(j)=0.
1115	fs(j)=w1
1116	endif
1117	endif
1118
1119	20 continue
1120
1121	* conditions aux limites pour les flux aux interfaces
1122
1123	fs(1)=0.
1124	ft(1)=0.
1125	fs(nz+1)=0.
1126	ft(nz+1)=-w1
1127
1128	* calcul des coefficients de l'equation discrete
1129
1130	do 21 j=nztop,nz
1131	as(j)=-dz(j)/dt
1132	bs(j)=-ft(j)
1133	cs(j)=ft(j+1)+fs(j)-dz(j)/dt
1134	ds(j)=-fs(j+1)
1135
1136	if ( cs(j).gt.0) goto 100
1137	21 continue
1138
1139	* cas explicite (mu=0) : calcul de la fonction c(z,r,t+1)
1140
1141	do 22 j=nztop,nz-1
1142
1143	if (j.eq.nztop) then
1144	dc=(cs(nztop)*c(nztop,k,li)+ds(nztop)
1145	& *c(nztop+1,k,li))/as(nztop)
1146	c(nztop,k,lf)=dc
1147	goto 22
1148	endif
1149
1150	dc=(bs(j)c(j-1,k,li)+cs(j)c(j,k,li)+ds(j)*c(j+1,k,li))
1151	s /as(j)
1152	c(j,k,lf)=dc
1153
1154
1155	22 continue
1156
1157	dc=(bs(nz)c(nz-1,k,li)+cs(nz)c(nz,k,li))/as(nz)
1158	c(nz,k,lf)=dc
1159
1160	if (nztop.ne.1) then
1161	do 32 j=1,nztop-1
1162	c(j,k,lf)=c(j,k,li)
1163	32 continue
1164	endif
1165
1166	goto 10
1167
1168	100 continue
1169
1170	* cas implicite (mu=1) : calcul de la fonction c(z,r,t+1)
1171
1172	do 101 j=nztop,nz
1173	asi(j)=ft(j)
1174	bsi(j)=-(ft(j+1)+fs(j)+dz(j)/dt)
1175	csi(j)=fs(j+1)
1176	dsi(j)=-dz(j)/dt*c(j,k,li)
1177	xsol(j)=0.
1178	101 continue
1179
1180	* inversion de la matrice tridiagonale
1181
1182	nb=nz-nztop+1
1183
1184	call dtridgl(nb,asi,bsi,csi,dsi,xsol)
1185
1186
1187	do 102 j=nztop,nz
1188	c(j,k,lf)=xsol(j)
1189	102 continue
1190
1191	if (nztop.ne.1) then
1192	do 110 j=1,nztop-1
1193	c(j,k,lf)=c(j,k,li)
1194	110 continue
1195	endif
1196
1197
1198
1199	10 continue
1200
1201	bilan14=0.
1202	do k=1,nrad
1203	do j=nztop,nz
1204	bilan14=bilan14+c(j,k,lf)dzb(j)
1205	& 4./3.pirf(k)*3.vrat_e**(k-imono)
1206	enddo
1207	enddo
1208
1209
1210	dice3=(bilan14-bilan4)*rhol
1211
1212
1213	return
1214
1215	end
1216

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

Download in other formats:

Original Format