Context Navigation

brume3D.F @ 1243

Last change on this file since 1243 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 30.6 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
[1056]	52	real dt,g
[175]	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
[1056]	89	* initialisation unique
	90	* --------------
	91
	92	c if (itime.eq.0) then
	93	c ITIME=1
	94	c endif
	95
	96	* initialisation
	97	* --------------
	98
[175]	99	call init
	100	call calcoag
	101
[1056]	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
[175]	116	do i=1,nz,1
	117	do j=1,nrad
	118	v1=vitesse(i,j,0)
[1056]	119	g=g0(rtit/(rtit+z(i)))*2
	120	taused(i,j)=rgpt(i)/(mn2g)/v1
[175]	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)
[474]	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
[175]	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.
[1056]	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
[175]	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
[1056]	769	rap=0.02
	770	c rap=0.191
[175]	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
[1056]	781	do i=1,nz
[175]	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)
[1056]	785	enddo
[175]	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