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

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)
	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