Context Navigation

brume3D.F @ 881

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

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

Download in other formats:

Original Format