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

Rev	Line
[175]	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)
[474]	115	c ho que c'est moche ! -> taused = RT/(Mn2g)1/vaer = H/v
[175]	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)
[474]	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
[175]	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