Context Navigation

brume3D.F @ 1643

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

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

Download in other formats:

Original Format