Context Navigation

snuages3D.F @ 1242

Last change on this file since 1242 was 1056, checked in by slebonnois, 11 years ago
SL: Titan runs ! see DOC/chantiers/commit_importants.log
File size: 46.6 KB

Rev	Line
[175]	1	subroutine snuages(ngrid,tab1,tab2,tab3,tab4,tab5,
	2	& x1,xnz,xnrad,ihor,fluxi,taused,precip)
	3
	4	!*
	5	!
	6	!1 c quantite de noyaux de la grille de rayon r a l'altitude z
	7	!2 c quantite de glace 1 de la grille de rayon r a l'altitude z
	8	!3 c quantite de glace 2 de la grille de rayon r a l'altitude z
	9	!4 c quantite de glace 3 de la grille de rayon r a l'altitude z
	10	!5 c quantite d'aerosols de la grille de rayon r a l'altitude z
	11	c
	12	c
	13	!x1 timestep.
	14
	15	c Notes : taused n'est actuellement pas calculé !
	16	c taused = dz/vitesse(z)
	17	c
	18	c
	19	--------------------------------------------------------------
	20	* *
	21	* ENTRE 0 ET 1000 KILOMETRES *
	22	* *
	23	* la dimension fractale est en tableau, attention au *
	24	* raccordement entre le regime moleculaire et le regime *
	25	* fluide *
	26	* *
	27	* Modele microphysique: Cabane et al.,1992 / *
	28	* Modele version fractale: Cabane et al.,1993 / *
	29	* *
	30	--------------------------------------------------------------
	31	* VERSION DU 2 JUIN 1993 --- AUT 1994 --- 11/04/96
	32	*
	33	* changer: altitude de production z0=/taux de production ctot=
	34	* : la charge/micron, ne
	35	* : df(h),rf...
	36	* raccordement aknc
	37	*
	38	* declaration des blocs communs
	39	*------------------------------
	40
	41	use dimphy
	42	IMPLICIT NONE
	43	#include "dimensions.h"
	44	#include "microtab.h"
	45	#include "varmuphy.h"
	46
	47	common/x2ctps/li,lf,dt
	48	common/x2con/c,c1,c2,c3,caer
	49	common/cldpart/rad,mmu
	50	common/x2frac/rfg,dfg
	51	common/x2effets/ xsaison
	52
	53
	54	* declaration des variables communes
	55	* ----------------------------------
	56
	57	integer xnz,xnrad,ngrid
	58	integer i,j,k,ihor,ibidon
	59	integer ihormx,imx
	60	integer li,lf
	61	real xsaison
	62	real dt
	63	real c(nz,nrad,2)
	64	real c1(nz,nrad,2),c2(nz,nrad,2)
	65	real c3(nz,nrad,2)
	66	real caer(nz,nrad,2)
	67	real rad(nz,nrad), mmu(nz,nrad)
	68	real rtemp(nz,nrad),mmutemp(nz,nrad)
	69	real knu,knu2,w
	70	real rfg(nz),dfg(nz,nrad)
	71	real ddt,dtau,vitesse2
	72	real vmax,vmin,rmax,rmin
	73	real fluxi(ngrid,nz,3)
	74	real taused(ngrid,nz)
	75	real precip(ngrid,5)
	76	real rmx
	77
	78
	79	* variables internes
	80	* ------------------
	81
	82	integer h,ti,itime,icoal
	83	real tab1(nz,nrad),tab2(nz,nrad),
	84	& tab3(nz,nrad),tab4(nz,nrad)
	85	real tab5(nz,nrad)
	86	real x1,xvolume,xmasse,xnoyaux
	87	real knd1,knd2,knd3,knd4,knd5,knd6
	88
	89	real dice1,dice2,dice3,dice4
	90
	91	data itime/0/
	92	data icoal/0/
	93	save itime,icoal
	94	save vmax,vmin,rmax,rmin
	95
	96
	97	* controles
	98	* ---------
	99
	100	if (nrad.ne.xnrad) stop 'nrad.ne.xnrad'
	101	if (nz.ne.xnz) stop 'nz.ne.xnz'
	102
	103	do i=1,nz
	104	do j=1,nrad
	105	c(i,j,1)=tab1(i,j)
	106	c(i,j,2)=0.0
	107	c1(i,j,1)=tab2(i,j)
	108	c1(i,j,2)=0.0
	109	c2(i,j,1)=tab3(i,j)
	110	c2(i,j,2)=0.0
	111	c3(i,j,1)=tab4(i,j)
	112	c3(i,j,2)=0.0
	113	caer(i,j,1)=tab5(i,j)
	114	caer(i,j,2)=0.0
	115	enddo
	116	enddo
	117
	118
	119	* initialisation
	120	* --------------
	121
	122	* Parametres physiques de Titan
	123
	124	if (itime.eq.0) then
	125	call init2
	126	print*,'init2 dans snuages.F'
	127	itime=1
	128	endif
	129
	130	aknc=2.92 !<--------Df=3
	131
	132	* Propriete des gouttes nuageuses/ calculee ici pour la sedimentation
	133	* r(i,j)=r (cst pour la colonne). Cela evite les accumulations dans
	134	* certains niveaux du a des differences de vitesse de sedimentation...
	135
	136	* ATTENTION, CES DEFINITIONS SONT INTERNES/ VOIR MUPHYS.F POUR DEFINITIONS
	137	* EXTERNES (UTILISEES DANS OPTCV ET OPTCI, PAR EX.)
	138
	139	do i=1,nz
	140	xnoyaux=0.
	141	xvolume=0.
	142	xmasse=0.
	143	do j=1,nrad
	144	dfg(i,j)=3.00
	145	rfg(j)=6.6e-8 !<--- arbitraire pour df=3
	146	xnoyaux=xnoyaux+c(i,j,1)
	147	xvolume=xvolume+c1(i,j,1)+c2(i,j,1)+c3(i,j,1)+
	148	& v_e(j)*c(i,j,1)
	149	xmasse =xmasse+
	150	& c1(i,j,1)*rhoi_ch4+
	151	& c2(i,j,1)*rhoi_c2h6+
	152	& c3(i,j,1)*rhoi_c2h2+
	153	& v_e(j)c(i,j,1)rhol
	154	enddo
	155
	156	do j=1,nrad
	157
	158	if (xnoyaux .le. 1.e-25 .or. xvolume.eq.0.) then ! utile ?
	159	dfg(i,j)=3.00
	160	rtemp(i,j) = 1.e-9
	161	mmutemp(i,j)= rhol !on prend la masse vol des aerosols
	162	else
	163	c Mais pourquoi se compliquer la vie alors que de toute facon on
	164	c reforce la masse volumique a son veritable calcul ^^
	165	c
	166	if ((c1(i,j,1)+c2(i,j,1)+c3(i,j,1))/xvolume.le.0.1) then ! glace / total
	167	dfg(i,j)=3.00
	168	rtemp(i,j) = ( (xvolume/xnoyaux)0.75/pi)*(1./3.)
	169	c mmutemp(i,j) = 1000.
	170	mmutemp(i,j)= xmasse/xvolume
	171	if(rtemp(i,j).gt.3.e-4) rtemp(i,j)=3.e-4
	172	else
	173	dfg(i,j)=3.00
	174	rtemp(i,j) = ( (xvolume/xnoyaux)0.75/pi)*(1./3.)
	175	mmutemp(i,j)= xmasse/xvolume
	176	if(rtemp(i,j).gt.3.e-4) rtemp(i,j)=3.e-4
	177	endif
	178	endif
	179	enddo
	180	enddo
	181
	182	do i=1,nz
	183	do j=1,nrad
	184	rad(i,j)=rtemp(i,1) ! mm valeur qqsoit j
	185	mmu(i,j)=mmutemp(i,1)
	186	enddo
	187	enddo
	188
	189
	190	* Coefficients de coagulation
	191
	192	if (icoal.eq.1) call calcoag2 ! 2 <-- 1
	193	vmin=1.e+6
	194	vmax=0.
	195	rmin=1.e+6
	196	rmax=0.
	197	ihormx=1
	198	imx=1
	199
	200
	201	* choix interne du temps d iteration
	202	* ----------------------------------
	203
	204	dt=x1
	205
	206	* iteration du modele sur le temps
	207	* ---------------------------------
	208
	209	li=1
	210	lf=2
	211
	212	* li=1 lf=2
	213
	214	call sedifn_fast(ihor,dice1,dice2,dice3,dice4) ! 1 <-- 1
	215	! call sedifn ! 2 <-- 1
	216
	217	li=3-li ! li devient 2
	218	lf=3-lf ! lf devient 1
	219
	220	if (icoal.eq.1) then
	221	* li=2 lf=1
	222	call coagul2(ihor) ! 1 <-- 2
	223
	224	li=3-li ! li devient 1
	225	lf=3-lf ! lf devient 2
	226
	227	endif
	228
	229	do i=1,nz
	230	vmin=vitesse2(i,1,1)
	231	do j=1,nrad
	232	fluxi(ihor,nz+1-i,1)=fluxi(ihor,nz+1-i,1)
	233	& -vminc1(i,j,li)rhoi_ch4
	234	fluxi(ihor,nz+1-i,2)=fluxi(ihor,nz+1-i,2)
	235	& -vminc2(i,j,li)rhoi_c2h6
	236	fluxi(ihor,nz+1-i,3)=fluxi(ihor,nz+1-i,3)
	237	& -vminc3(i,j,li)rhoi_c2h2
	238	enddo
	239	enddo
	240
	241	c En theorie, les diceX sont NEGATIF (en sedimentant on ne fait que perdre de la glace)
[474]	242	c Les precipitations sont comptees positivement. (ET ON NE PREND QUE DES VALEURS POSITIVES !)
[175]	243
[814]	244	precip(ihor,1)=AMAX1(-dice1/rhoi_ch4,0.) ! m3/m2 = m :)
	245	precip(ihor,2)=AMAX1(-dice2/rhoi_c2h6,0.)
	246	precip(ihor,3)=AMAX1(-dice3/rhoi_c2h2,0.)
	247	precip(ihor,4)=AMAX1(-dice4/rhol,0.)
[474]	248
[175]	249	do i=1,nz
	250	do j=1,nrad
	251	tab1(i,j)=c(i,j,li) ! li=1
	252	tab2(i,j)=c1(i,j,li) ! li=1
	253	tab3(i,j)=c2(i,j,li) ! li=1
	254	tab4(i,j)=c3(i,j,li) ! li=1
	255	enddo
	256	enddo
	257
	258	return
	259
	260	end
	261
	262
	263
	264	*______________________________________________________________________
	265
	266	real function lambda2(j,indic)
	267	*
	268	------------------------------------------------------------------
	269	* fonction calculant le libre parcours moyen des molecules *
	270	* atmospheriques( rayon =ra) se trouvant dans la couche no j. *
	271	* pour indic=0 ...... la particule se trouve a la frontiere entre*
	272	* les couches j et j-1 *
	273	* pour indic=1 ...... la particule se trouve au milieu de la *
	274	* la couche j *
	275	------------------------------------------------------------------
	276	*
	277	* declaration des blocs communs
	278	*------------------------------
	279	use dimphy
	280	IMPLICIT NONE
	281	#include "dimensions.h"
	282	#include "microtab.h"
	283	#include "varmuphy.h"
	284
	285	common/x2frac/rfg,dfg
	286
	287	* declaration des variables communes
	288	* ----------------------------------
	289
	290	real rfg(nz),dfg(nz,nrad)
	291
	292	* declaration des variables internes
	293	* ----------------------------------
	294
	295	integer indic,j
	296	real pp,ra
	297
	298	ra=1.75e-10
	299
	300	* traitement
	301	* ----------
	302
	303	if (indic.eq.0) then
	304	pp=pb(j)
	305	else
	306	if (indic.ne.1) then
	307	print*,'erreur argument fonction lambda'
	308	stop
	309	endif
	310	pp=p(j)
	311	endif
	312
	313	lambda2=kbzt(j)/(4sqrt(2.)pi(ra*2)pp)
	314	end
	315
	316	*******************************************************************************
	317
	318	real function knu2(j,k,indic)
	319	*
	320	--------------------------------------------------------------
	321	* fonction calculant le nombre de knudsen d'une particule *
	322	* d'aerosol de rayon rad(k) se trouvant dans la couche no j *
	323	* indic ...... idem function lambda *
	324	--------------------------------------------------------------
	325	*
	326	* declaration des blocs communs
	327	*------------------------------
	328	use dimphy
	329	IMPLICIT NONE
	330	#include "dimensions.h"
	331	#include "microtab.h"
	332	#include "varmuphy.h"
	333
	334	common/cldpart/rad,mmu
	335	common/x2frac/rfg,dfg
	336
	337
	338	* declaration des variables communes
	339	* ----------------------------------
	340
	341	real rad(nz,nrad),mmu(nz,nrad)
	342	real rfg(nz),dfg(nz,nrad)
	343
	344
	345	* declaration des variables internes
	346	* ----------------------------------
	347
	348	integer indic,j,k
	349	real lambda2,rfk
	350
	351	* traitement
	352	* ----------
	353
	354	if (indic.ne.0 .and.indic.ne.1) then
	355	print*,'erreur argument fonction knu'
	356	stop
	357	endif
	358
	359
	360	rfk=(rad(j,k)*(3./dfg(j,k)))((rfg(k))**(1.-3./dfg(j,k)))
	361	knu2=lambda2(j,indic)/rfk
	362
	363	end
	364
	365	*****************************************************************************
	366
	367	real function nud2(j,indic)
	368	*
	369	--------------------------------------------------------------
	370	* fonction calculant la viscosite dynamique (en USI) de l air *
	371	* d apres la formule de Sutherlant a l altitude j *
	372	* indic ......... idem fonction lambda *
	373	--------------------------------------------------------------
	374	*
	375	use dimphy
	376	IMPLICIT NONE
	377	#include "dimensions.h"
	378	#include "microtab.h"
	379	#include "varmuphy.h"
	380	integer indic,j
	381	real nud0,c,tt
	382	real rfg(nz),dfg(nz,nrad)
	383
	384	common/x2frac/rfg,dfg
	385
	386	*
	387	nud0=1.74e-5
	388	c=109.
	389
	390	if(indic.ne.0.and.indic.ne.1) then
	391	print*,'erreur argument fonction nud'
	392	stop
	393	endif
	394
	395	if (indic.eq.0) tt=tb(j)
	396	if (indic.eq.1) tt=t(j)
	397	nud2=nud0sqrt(tt/293)(1+c/293)/(1+c/tt)
	398	end
	399
	400	****************************************************************************
	401
	402	real function vitesse2(j,k,indic)
	403	*
	404	-----------------------------------------------------------------
	405	* fonction calculant la vitesse de chute d une particule de rayon*
	406	* k se trouvant a l altitude j suivant la valeur du nombre de *
	407	* Knudsen *
	408	* indic ....... idem function lambda *
	409	-----------------------------------------------------------------
	410	*
	411
	412	* declaration des blocs communs
	413	*------------------------------
	414	use dimphy
	415	IMPLICIT NONE
	416	#include "dimensions.h"
	417	#include "microtab.h"
	418	#include "varmuphy.h"
	419	common/cldpart/rad,mmu
	420	common/x2frac/rfg,dfg
	421	common/x2ctps/li,lf,dt
	422
	423	* declaration des variables communes
	424	* ----------------------------------
	425
	426	real rad(nz,nrad),mmu(nz,nrad)
	427	real rfg(nz),dfg(nz,nrad)
	428
	429
	430	* declaration des variables internes
	431	* ----------------------------------
	432
	433	integer indic,j,k,li,lf
	434	real w,g,m,a0,zz,knu2,nud2,knud,tt,rhoh
	435	real vlimite,akncx,rbis,rfk
	436	real dt
	437
	438	* traitement
	439	* ----------
	440
	441	if (indic.ne.0.and.indic.ne.1) then
	442	print*,'erreur argument fonction vitesse'
	443	stop
	444	endif
	445
	446	if(indic.eq.0) then
	447	zz=z(j)+dz(j)/2.
	448	tt=tb(j)
	449	rhoh=rhob(j)
	450	endif
	451	if(indic.eq.1) then
	452	zz=z(j)
	453	tt=t(j)
	454	rhoh=rho(j)
	455	endif
	456
	457	g=g0(rtit/(rtit+zz))*2
	458	a0=0.74
	459	m=(ach4(j)mch4+aar(j)mar+an2(j)*mn2)/nav
	460	knud=knu2(j,k,indic)
	461
	462	c akncx=aknc
	463	c if(df(k).gt.2.5) akncx=2.7
	464
	465	c if(knud.ge.akncx) then
	466	c rbis=(rad(j,k)*(3.-6./dfg(j,k)))((rfg(k))**(-2.+6./dfg(j,k)))
	467	c w=a0grbismmu(j,k)/(rhohsqrt(8kbztt/(pi*m)))
	468	c endif
	469
	470	rfk=(rad(j,k)*(3./dfg(j,k)))((rfg(k))**(1.-3./dfg(j,k)))
	471	w=2./9.rfk(dfg(j,k)-1.)rfg(k)*(3.-dfg(j,k))g*mmu(j,k)
	472	& /nud2(j,indic)
	473
	474	w=w(1+1.2517knud+0.4knudexp(-1.1/knud))
	475
	476	w=w!*3. ! on tient compte de la largeur de distribution... a affiner
	477	vitesse2=w
	478
	479
	480	end
	481	***********************************************************************
	482	real function kd2(h)
	483	*
	484	--------------------------------------------------------------------
	485	* cette fonction calcule le coefficient du terme de eddy diffusion *
	486	* a l altitude j *
	487	--------------------------------------------------------------------
	488	*
	489	use dimphy
	490	IMPLICIT NONE
	491	#include "dimensions.h"
	492	#include "microtab.h"
	493	#include "varmuphy.h"
	494
	495	common/x2frac/rfg,dfg
	496
	497	real zbx
	498	real rfg(nz),dfg(nz,nrad)
	499
	500	integer h
	501
	502	zbx=z(h)+dz(h)/2.
	503	if(zbx.le.42000.) then
	504	c kd2=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
	505	kd2=4.
	506	else
	507	kd2=1.64e+12(pb(h)/(kbztb(h)))**(-1./2.)
	508	kd2=0.0*kd2
	509	endif
	510	kd2=00.
	511
	512
	513	return
	514	end
	515
	516
	517	*____________________________________________________________________________
	518
	519	subroutine init2
	520	*
	521	--------------------------------------------------------------------
	522	* cette routine effectue : *
	523	* 1) interpolation a partir des donnees initiales des *
	524	* valeurs de p,t,rho,ach4,aar,an2 sur la grille *
	525	* 2) initialisation des constantes (common/x2phys/) *
	526	* 3) initialisation des variables temporelles (common *
	527	* /temps/) *
	528	* 4) definition des grilles en rayon et verticale *
	529	* 5) initialisation de c(z,r,t) avec les donnees du *
	530	* fichier unit=1 *
	531	* *
	532	* les donnees sont des valeurs caracterisques de l atmosphere de *
	533	* TITAN ( voir Lelouch and co ) *
	534	--------------------------------------------------------------------
	535
	536	* declaration des blocs communs
	537	*------------------------------
	538	use dimphy
	539	IMPLICIT NONE
	540	#include "dimensions.h"
	541	#include "microtab.h"
	542	#include "varmuphy.h"
	543	common/x2ctps/li,lf,dt
	544	common/cldpart/rad,mmu
	545
	546	* declaration des variables communes
	547	* ----------------------------------
	548
	549	integer li,lf
	550	real dt
	551	real rad(nz,nrad), mmu(nz,nrad)
	552
	553
	554	* declaration des variables internes
	555	* ----------------------------------
	556	integer nzd,i,ii
	557	parameter (nzd=254)
	558	integer limsup,liminf,j1,j2
	559	real zd(nzd),ach4d(nzd),rap
	560	real m
	561
	562
	563	* initialisation des variables temporelles
	564	* ----------------------------------------
	565
	566	li=1
	567	lf=2
	568
	569
	570	* interpolation de xch4,xar et xn2 sur la grille
	571	* ----------------------------------------------
	572
	573	* donnees initiales (Lellouch et al,87)
	574	* -------------------------------------
	575
	576	c print,'***** init'
	577	do 1 i=1,168
	578	zd(i)=(1000.-5(i-1))1000.
	579	1 continue
	580	do 2 i=1,78
	581	zd(168+i)=(160.-2(i-1))1000.
	582	2 continue
	583	do 3 i=1,4
	584	zd(246+i)=(5.-(i-1))*1000.
	585	3 continue
	586	do 4 i=1,4
	587	zd(250+i)=(1.5-(i-1)0.5)1000.
	588	4 continue
	589
	590	data (ach4d(i),i=1,168)/168*1.5e-2/
	591	data (ach4d(i),i=169,254)/63*1.5e-2,1.6e-2,1.8e-2,1.8e-2,
	592	& 1.9e-2,2.e-2,2.1e-2,2.3e-2,2.5e-2,2.8e-2,3.1e-2,3.6e-2,
	593	& 4.1e-2,4.7e-2,5.7e-2,6.7e-2,7.5e-2,7*8.e-2/
	594
	595	liminf=0
	596	limsup=0
	597
	598	* interpolation des taux de melange de ch4,ar,n2
	599	*-----------------------------------------------
	600
	601	do 20 j1=1,nz
	602	do 21 j2=1,nzd
	603	if( zd(j2).le.z(j1)) goto 22
	604	21 continue
	605	22 liminf=j2
	606	if (zd(liminf).eq.z(j1) )then
	607	ach4(j1)=ach4d(liminf)
	608	goto 20
	609	endif
	610	if (j2.ne.1) then
	611	limsup=j2-1
	612	else
	613	limsup=j2
	614	endif
	615	if (limsup.eq.liminf) then
	616	ach4(j1)=ach4(limsup)
	617	else
	618	ach4(j1)=ach4d(liminf)-(ach4d(limsup)-ach4d(liminf))/
	619	s (zd(limsup)-zd(liminf))*(zd(liminf)-z(j1))
	620	endif
	621	20 continue
	622
	623	* rap= aar/an2 cst sur l altitude
	624
	625	rap=0.191
	626	do 23 i=1,nz
	627	an2(i)=(1.-ach4(i))/(1.+rap)
	628	aar(i)=rap*an2(i)
	629	23 continue
	630
	631	do 24 i=1,nz
	632	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
	633	rho(i)=p(i)m/(rgpt(i))
	634	24 continue
	635
	636	do 34 i=1,nz
	637	m=ach4(i)mch4+an2(i)mn2+aar(i)*mar
	638	rhob(i)=pb(i)m/(rgptb(i))
	639	c print*,pb(i),m,rgp,tb(i),rhob(i),rho(i)
	640	34 continue
	641
	642	* fin d interpolation des taux de melange
	643	*----------------------------------------
	644
	645	c print,'*** fin init'
	646	540 continue
	647	return
	648
	649	500 print*,'erreur lecture initialisation de c...erreur=',ii
	650	stop
	651
	652	end
	653
	654	*__________________________________________________________________________
	655
	656	subroutine sedifn
	657	*
	658	------------------------------------------------------------------
	659	* cette routine calcule l evolution de la fonction de distribution*
	660	* c(z,r,t) pour les phenomenes de sedimentation et de diffusion *
	661	------------------------------------------------------------------
	662	*
	663	*
	664	* declaration des blocs communs
	665	*------------------------------
	666	use dimphy
	667	IMPLICIT NONE
	668	#include "dimensions.h"
	669	#include "microtab.h"
	670	#include "varmuphy.h"
	671
	672	common/x2ctps/li,lf,dt
	673	common/x2con/ctmp,c1,c2,c3,caer
	674	common/cldpart/rad,mmu
	675	common/x2frac/rfg,dfg
	676
	677
	678	* declaration des variables communes
	679	* ----------------------------------
	680
	681	integer li,lf
	682	real dt
	683	real c(nz,nrad,2)
	684	real c1(nz,nrad,2)
	685	real c2(nz,nrad,2)
	686	real c3(nz,nrad,2)
	687	real caer(nz,nrad,2)
	688	real ctmp(nz,nrad,2)
	689	real rad(nz,nrad),mmu(nz,nrad)
	690	real rfg(nz),dfg(nz,nrad)
	691
	692
	693	* declaration des variables internes
	694	* ----------------------------------
	695
	696	real w,w1,dzbX,dc
	697	integer i,j,k,ii,nb
	698	double precision sigma,theta,hc,l,rap,cmp,wp
	699	double precision fs(nz+1),ft(nz+1)
	700	real as(nz),bs(nz),cs(nz),ds(nz)
	701	double precision asi(nztop:nz),bsi(nztop:nz),
	702	& csi(nztop:nz)
	703	double precision dsi(nztop:nz),xsol(nztop:nz)
	704	real vitesse2,kd2
	705
	706	external dtridgl
	707
	708	* resolution
	709	*------------
	710
	711	do 5 ii=1,4
	712	do k=1,nrad
	713	do j=nztop,nz
	714	if( ii.eq.1 ) c(j,k,li)=ctmp(j,k,li)
	715	if( ii.eq.2 ) c(j,k,li)=c1(j,k,li)
	716	if( ii.eq.3 ) c(j,k,li)=c2(j,k,li)
	717	if( ii.eq.4 ) c(j,k,li)=c3(j,k,li)
	718	enddo
	719	enddo
	720
	721	do 10 k=1,nrad
	722	do 20 j=nztop,nz
	723	if (j.eq.1) goto 20
	724	* calcul de la vitesse corrigee
	725
	726	dzbX=(dz(j)+dz(j-1))/2.
	727	w= -1*vitesse2(j,k,0)
	728
	729	if (kd2(j).ne.0.) then
	730	theta=0.5(wdzbX/kd2(j)+log(rho(j-1)/rho(j)))
	731	if (theta.ne.0) then
	732	sigma=1./dtanh(theta)-1./theta
	733	else
	734	sigma=1.
	735	endif
	736	else
	737	sigma=1.
	738	endif
	739
	740	if(c(j,k,li).eq.0.) then
	741	rap=10.
	742	else
	743	rap=c(j-1,k,li)/c(j,k,li)
	744	if( rap.gt.10.) rap=10.
	745	if( rap.lt.0.1) rap=0.1
	746	endif
	747
	748	if (rap.gt.0.9 .and. rap.lt.1.1) then
	749	w1=w
	750	else
	751	if(w.ne.0) then
	752	hc=dzbX/dlog(rap)
	753	l=dzbX/(wdt)(dexp(-w*dt/hc)-1.)/(1.-rap)
	754	wp=w*1.d0
	755	cmp=dlog(-wp)+abs(sigma)*dlog(l)
	756	if (cmp.gt.38) then
	757	goto 20
	758	endif
	759	w1=-dexp(cmp)
	760	else
	761	w1=0.
	762	endif
	763	endif
	764
	765	* calcul des flux aux interfaces
	766
	767
	768	if (kd2(j).ne.0.) then
	769	if (theta.ne.0.) then
	770	ft(j)=(w1+log(rho(j-1)/rho(j))*kd2(j)/dzbX)/
	771	& (dexp(2.*theta)-1.)
	772	fs(j)=ft(j)dexp(2.theta)
	773	else
	774	ft(j)=kd2(j)/dzbX
	775	fs(j)=kd2(j)/dzbX
	776	endif
	777	else
	778	if (w1.lt.0.)then
	779	ft(j)=-w1
	780	fs(j)=0.
	781	else
	782	ft(j)=0.
	783	fs(j)=w1
	784	endif
	785	endif
	786
	787	20 continue
	788
	789	* conditions aux limites pour les flux aux interfaces
	790
	791	fs(1)=0.
	792	ft(1)=0.
	793	fs(nz+1)=0.
	794	ft(nz+1)=-w1
	795
	796	* calcul des coefficients de l equation discrete
	797
	798	do 21 j=nztop,nz
	799	as(j)=-dz(j)/dt
	800	bs(j)=-ft(j)
	801	cs(j)=ft(j+1)+fs(j)-dz(j)/dt
	802	ds(j)=-fs(j+1)
	803	if ( cs(j).gt.0) goto 100
	804	21 continue
	805
	806	* cas explicite (mu=0) : calcul de la fonction c(z,r,t+1)
	807
	808	do 22 j=nztop,nz-1
	809
	810	if (j.eq.nztop) then
	811	dc=(cs(nztop)*c(nztop,k,li)+ds(nztop)
	812	& *c(nztop+1,k,li))/as(nztop)
	813	c(nztop,k,lf)=dc
	814	goto 22
	815	endif
	816
	817	dc=(bs(j)c(j-1,k,li)+cs(j)c(j,k,li)+ds(j)*c(j+1,k,li))
	818	s /as(j)
	819	c(j,k,lf)=dc
	820
	821	22 continue
	822
	823	dc=(bs(nz)c(nz-1,k,li)+cs(nz)c(nz,k,li))/as(nz)
	824	c(nz,k,lf)=dc
	825
	826	if (nztop.ne.1) then
	827	do 32 j=1,nztop-1
	828	c(j,k,lf)=c(j,k,li)
	829	32 continue
	830	endif
	831
	832	goto 10
	833
	834	100 continue
	835
	836	* cas implicite (mu=1) : calcul de la fonction c(z,r,t+1)
	837
	838	do 101 j=nztop,nz
	839	asi(j)=ft(j)
	840	bsi(j)=-(ft(j+1)+fs(j)+dz(j)/dt)
	841	csi(j)=fs(j+1)
	842	dsi(j)=-dz(j)/dt*c(j,k,li)
	843	101 continue
	844
	845	* inversion de la matrice tridiagonale
	846
	847	nb=nz-nztop+1
	848
	849	call dtridgl(nb,asi,bsi,csi,dsi,xsol)
	850
	851	do 102 j=nztop,nz
	852	c(j,k,lf)=xsol(j)
	853	102 continue
	854
	855	if (nztop.ne.1) then
	856	do 110 j=1,nztop-1
	857	c(j,k,lf)=c(j,k,li)
	858	110 continue
	859	endif
	860
	861
	862
	863	10 continue
	864
	865	do k=1,nrad
	866	do j=nztop,nz
	867	if( ii.eq.1 ) ctmp(j,k,lf)=c(j,k,lf)
	868	if( ii.eq.2 ) c1(j,k,lf) =c(j,k,lf)
	869	if( ii.eq.3 ) c2(j,k,lf) =c(j,k,lf)
	870	if( ii.eq.4 ) c3(j,k,lf) =c(j,k,lf)
	871	enddo
	872	enddo
	873
	874	5 continue
	875
	876	return
	877
	878	end
	879
	880	*__________________________________________________________________________
	881
	882	subroutine sedifn_fast(ihor,dice1,dice2,dice3,dice4)
	883	*
	884	------------------------------------------------------------------
	885	* cette routine calcule l evolution de la fonction de distribution *
	886	* c(z,r,t) pour les phenomenes de sedimentation {pas de diffusion} *
	887	* dice1 = delta glace CH4 *
	888	* dice2 = delta glace C2H6 *
	889	* dice3 = delta glace C2H2 *
	890	* dice4 = delta Volume noyaux *
	891	------------------------------------------------------------------
	892	*
	893	* declaration des blocs communs
	894	*------------------------------
	895	use dimphy
	896	IMPLICIT NONE
	897	#include "dimensions.h"
	898	#include "microtab.h"
	899	#include "varmuphy.h"
	900
	901	common/x2ctps/li,lf,dt
	902	common/x2con/c,c1,c2,c3,caer
	903	common/cldpart/rad,mmu
	904	common/x2frac/rfg,dfg
	905
	906
	907	* declaration des variables communes
	908	* ----------------------------------
	909
	910	integer li,lf,i,j,k,ihor
	911	integer jinf,jsup,jj,iiter
	912	real dt
	913	real c(nz,nrad,2)
	914	real c1(nz,nrad,2)
	915	real c2(nz,nrad,2)
	916	real c3(nz,nrad,2)
	917	real caer(nz,nrad,2)
	918	real ci(nz,nrad,2)
	919	real ci1(nz,nrad,2)
	920	real ci2(nz,nrad,2)
	921	real ci3(nz,nrad,2)
	922	real rad(nz,nrad),mmu(nz,nrad)
	923	real rfg(nz),dfg(nz,nrad)
	924	real puit(nz)
[1056]	925	c ------ echange est cree sur la taille maxi mais n'est utilisee
	926	c que sur la dim geree par le proc (klon ou jjm+1)
[175]	927	integer ngrid
[1056]	928	parameter (ngrid=(jjm-1)*iim+2) ! = taille maximum
[175]	929	real echange(nz,nz,ngrid)
[1056]	930	c pas genial mais vu que c est tres local, pas de soucis a priori en parallele.
[175]	931	real bilan1,bilan2,bilan3,bilan4,bilan5
	932	real bilan11,bilan12,bilan13,bilan14,bilan15
	933	real compte,compte2,xepl
	934	real dice1,dice2,dice3,dice4
	935
	936
	937
	938	* declaration des variables internes
	939	* ----------------------------------
	940
	941	real vitesse2,kd2
	942	real w,ws,wi,zs,zi,alpha,v0,deltazs,deltazi
	943	real zni,znip1,xcnt,ichx,arg1,arg2,xft,xf
	944	real argexp,explim
	945
	946	save echange
	947
	948
	949	* resolution
	950	*------------
	951
	952
	953	bilan1=0.
	954	bilan2=0.
	955	bilan3=0.
	956	bilan4=0.
	957	bilan5=0.
	958	do k=1,nrad
	959	do j=nztop,nz
	960	ci(j,k,li)= c(j,k,li)*dzb(j) ! li
	961	ci1(j,k,li)=c1(j,k,li)*dzb(j)
	962	ci2(j,k,li)=c2(j,k,li)*dzb(j)
	963	ci3(j,k,li)=c3(j,k,li)*dzb(j)
	964	bilan5=bilan5+ci(j,k,li)
	965	bilan1=bilan1+ci1(j,k,li)
	966	bilan2=bilan2+ci2(j,k,li)
	967	bilan3=bilan3+ci3(j,k,li)
	968	bilan4=bilan4+
	969	& ci(j,k,li)4./3.pirf(k)3.vrat_e**(k-imono)
	970
	971	ci(j,k,lf)= 0. ! lf
	972	ci1(j,k,lf)=0.
	973	ci2(j,k,lf)=0.
	974	ci3(j,k,lf)=0.
	975	enddo
	976	enddo
	977
	978	* calcul de la matrice d echange
	979	*----------------------------------------------------------------
	980
	981	do j=nztop,nz
	982	do i=nztop,nz
	983	echange(i,j,ihor)=0.
	984	enddo
	985	enddo
	986
	987	do 20 i=nztop,nz
	988	puit(i)=0.
	989	do 30 k=1,1
	990	xcnt=0.
	991	ICHX=1 ! extrapolation 0: lineaire 1: exponentielle
	992	IF(ICHX.eq.0 .or. ICHX.eq.2) THEN
	993	ws=vitesse2(i,k,0)
	994	if(i.lt.nz) wi=vitesse2(i+1,k,0)
	995	if(i.eq.nz) wi=vitesse2(i ,k,1)
	996	w=(ws+wi)/2.
	997	zni =zb(i)-w*dt
	998	znip1=zb(i)-dzb(i)-w*dt
	999	ENDIF
	1000
	1001	explim=30.
	1002
	1003	IF(ICHX.eq.1 .or. ICHX.eq.2) THEN
	1004	ws=vitesse2(i,k,0)
	1005	zs=zb(i)
	1006	wi=vitesse2(i,k,1)
	1007	zi=z(i)
	1008
	1009	c if(wi.eq.ws) wi=ws/1.001 ! sinon ca plante !
	1010	if(abs(wi-ws)/wi .le. 1.e-3) wi=ws/1.001 ! sinon ca plante !
	1011
	1012	if(wi.ne.0.) alpha= alog(ws/wi) /(zs-zi)
	1013	argexp=alpha*zs
	1014	if(argexp.lt.-explim) argexp=-explim
	1015	if(argexp.gt. explim) argexp=+explim
	1016	v0 = ws/exp(argexp)
	1017
	1018	argexp=alpha*zb(i)
	1019	if(argexp.lt.-explim) argexp=-explim
	1020	if(argexp.gt. explim) argexp=+explim
	1021	arg1=1.+v0alphaexp(argexp)*dt
	1022
	1023	argexp=alpha*(zb(i)-dzb(i))
	1024	if(argexp.lt.-explim) argexp=-explim
	1025	if(argexp.gt. explim) argexp=+explim
	1026	arg2=1.+v0alphaexp(argexp)*dt
	1027
	1028	iiter=0
	1029	101 continue
	1030
	1031	if (iiter.le.25) then
	1032	if(arg1.le.0. .or. arg2.le.0.) then
	1033	print*,ihor,i,iiter, 'ajustement vitesse',arg1,arg2
	1034	print*,ws,wi, ' w1 w2 anc valeurs'
	1035	print*,alpha, ' alpha anc valeurs'
	1036	print*,rad(i,k), 'r(j,k)'
	1037	print*,mmu(i,k), ' mmu(j,k)'
	1038	print*,t(i), ' t(j,k)'
	1039	print*,arg1,arg2, ' arg1 arg2 anc valeurs'
	1040	arg2=1+(arg2-1.)/2.
	1041	arg1=1+(arg1-1.)/2.
	1042	iiter=iiter+1
	1043	print*,arg1,arg2, ' arg1 arg2 nle valeurs'
	1044	goto 101
	1045	endif
	1046	else
	1047	stop
	1048	endif
	1049
	1050	if(arg1.lt.0.) print*,'arg1:',arg1
	1051	if(arg2.lt.0.) print*,'arg2:',arg2
	1052
	1053	deltazs=-alog(arg1)/alpha
	1054	deltazi=-alog(arg2)/alpha
	1055
	1056	zni =zb(i)+deltazs
	1057	znip1=zb(i)-dzb(i)+deltazi
	1058
	1059	ENDIF
	1060
	1061	if(zni.ne.znip1) xft=zni/(zni-znip1)
	1062	if(zni.eq.znip1 .and. zni.le.0.) xft=0.
	1063	if(zni.eq.znip1 .and. zni.gt.0.) then
	1064	xft=0.
	1065	print*,'zni..znip1', zni,znip1
	1066	endif
	1067
	1068	* Si des aerosols touchent le sol (zni < 0 ) alors on fixe
	1069	* le niveau a 0, et on elimine les aerosols correspondants
	1070	*-----------------------------------------------------------
	1071
	1072	if(znip1 .lt. 0.) znip1=0.
	1073	if(zni .lt. 0.) zni =0.
	1074
	1075	if(zni.le.0. .and. znip1 .le. 0.) then
	1076	c print*,'voie 1 / disparition complete'
	1077	xft=0.
	1078	xf=1.
	1079	xcnt=xcnt+xf
	1080	puit(i)=puit(i)+xf
	1081	endif
	1082
	1083	if(zni.gt.0. .and. znip1 .le. 0.) then
	1084	c print*,'voie 2 / disparitipon partielle'
	1085	xf=(1-xft)
	1086	xcnt=xcnt+xf
	1087	puit(i)=puit(i)+xf
	1088	endif
	1089
	1090	if(zni.gt.0. .and. znip1 .gt. 0.) then
	1091	c print*,'voie 3 / pas de disparition'
	1092	xft=1.
	1093	xf=(1-xft)
	1094	xcnt=xcnt+xf
	1095	puit(i)=puit(i)+xf
	1096	endif
	1097
	1098	jsup=nz+1
	1099	jinf=nz+1
	1100	do j=nztop,nz
	1101	if(zni.le.zb(j) .and. zni.ge.zb(j)-dzb(j)) jsup=j
	1102	if(znip1.le.zb(j).and. znip1.ge.zb(j)-dzb(j)) jinf=j
	1103	enddo
	1104	if(zni .ge. 0. .and. zni .lt. 1.e-3) jsup=nz
	1105	if(znip1 .ge. 0. .and. znip1 .lt. 1.e-3) jinf=nz
	1106
	1107
	1108	* Volume inclu dans un seul niveau
	1109	*----------------------------------
	1110
	1111	if (jsup .eq. jinf. and. jsup.ge. nz+1) then
	1112	xcnt=xcnt+1.
	1113	print*,'cas impossible'
	1114	print*,'alpha= ',alpha
	1115	print*,'ws wi ',ws,wi
	1116	print*,'deltazs deltazi ',deltazs,deltazi
	1117	print*,' r(i,k) mmu(i,k)', rad(i,k),mmu(i,k)
	1118	print*,' t(j,k)',t(i)
	1119	print*,zni,znip1,jsup,jinf
	1120	print*,'STOP'
	1121	STOP
	1122	endif
	1123
	1124	if (jsup .eq. jinf. and. jsup.le. nz) then
	1125	xf=1.
	1126	xcnt=xcnt+xft*xf
	1127	if(jinf.le.nz) then
	1128	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
	1129	endif
	1130	endif
	1131
	1132	* Volume a cheval sur 2 niveaux
	1133	*--------------------------------
	1134
	1135	if (jinf .eq. jsup+1) then
	1136	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
	1137	xcnt=xcnt+xf*xft
	1138	if(jsup.le.nz) then
	1139	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
	1140	endif
	1141	xf=(zb(jinf)-znip1)/(zni-znip1)
	1142	xcnt=xcnt+xf*xft
	1143	if(jinf.le.nz) then
	1144	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
	1145	endif
	1146	endif
	1147
	1148	* Volume a cheval sur 3 ou plus de niveaux
	1149	*------------------------------------------
	1150
	1151	if (jinf .gt. jsup+1) then
	1152
	1153	c print*,' voie C / dans N cases'
	1154	xf=(zni-zb(jsup)+dzb(jsup))/(zni-znip1)
	1155	xcnt=xcnt+xf*xft
	1156	if(jsup.le.nz) then
	1157	echange(jsup,i,ihor)=echange(jsup,i,ihor)+xft*xf
	1158	endif
	1159
	1160	xf=(zb(jinf)-znip1)/(zni-znip1)
	1161	xcnt=xcnt+xf*xft
	1162	if(jinf.le.nz) then
	1163	echange(jinf,i,ihor)=echange(jinf,i,ihor)+xft*xf
	1164	endif
	1165
	1166	do jj=jsup+1,jinf-1
	1167	xf=(zb(jj)-zb(jj+1))/(zni-znip1)
	1168	xcnt=xcnt+xf*xft
	1169	if(jj.le.nz) then
	1170	echange(jj,i,ihor)=echange(jj,i,ihor)+xft*xf
	1171	endif
	1172	enddo
	1173
	1174	endif
	1175
	1176
	1177	* et sur les rayons...
	1178	*---------------------
	1179	30 continue
	1180
	1181	* fin de la grande boucle sur les altitudes...
	1182	*----------------------------------------------
	1183	20 continue
	1184
	1185	* Calcul etat final Cfinal = Echange*initial
	1186	*----------------------------------------------
	1187
	1188	compte=0.
	1189	compte2=0.
	1190	do j=1,nz
	1191	xepl=0.
	1192	do jj=1,nz
	1193	xepl=xepl+echange(jj,j,ihor)
	1194	compte=compte+echange(jj,j,ihor)
	1195	compte2=compte2+echange(jj,j,ihor)
	1196	enddo
	1197	compte2=compte2+puit(j)
	1198	enddo
	1199
	1200	if(abs(compte2-nz) .gt. 1.e-4)
	1201	& print*,'Matrice calculee#',ihor,'tx d expl (/55):',
	1202	& compte,compte2
	1203
	1204
	1205	* Fin du calcul de la matrice d*echange
	1206	*----------------------------------------------
	1207
	1208	do j=nztop,nz
	1209	do k=1,nrad
	1210
	1211	do jj=nztop,nz
	1212	ci(j,k,lf)=ci(j,k,lf)+
	1213	& echange(j,jj,ihor)*ci(jj,k,li)
	1214	ci1(j,k,lf)=ci1(j,k,lf)+
	1215	& echange(j,jj,ihor)*ci1(jj,k,li)
	1216	ci2(j,k,lf)=ci2(j,k,lf)+
	1217	& echange(j,jj,ihor)*ci2(jj,k,li)
	1218	ci3(j,k,lf)=ci3(j,k,lf)+
	1219	& echange(j,jj,ihor)*ci3(jj,k,li)
	1220	enddo
	1221	enddo
	1222	enddo
	1223
	1224
	1225	* Controles et affichage Bilan
	1226	*----------------------------------------------
	1227
	1228	bilan11=0.
	1229	bilan12=0.
	1230	bilan13=0.
	1231	bilan14=0.
	1232	bilan15=0.
	1233	do k=1,nrad
	1234	do j=nztop,nz
	1235	c(j,k,lf) =ci(j,k,lf)/dzb(j)
	1236	c1(j,k,lf)=ci1(j,k,lf)/dzb(j)
	1237	c2(j,k,lf)=ci2(j,k,lf)/dzb(j)
	1238	c3(j,k,lf)=ci3(j,k,lf)/dzb(j)
	1239	bilan15=bilan15+ ci(j,k,lf)
	1240	bilan11=bilan11+ci1(j,k,lf)
	1241	bilan12=bilan12+ci2(j,k,lf)
	1242	bilan13=bilan13+ci3(j,k,lf)
	1243	bilan14=bilan14+
	1244	& ci(j,k,lf)4./3.pirf(k)3.vrat_e**(k-imono)
	1245	enddo
	1246	enddo
	1247
	1248	c print*,'sedifn_fast Bilans:'
	1249	c & ,bilan11,bilan12,bilan13
	1250	c print*,'Bilan1:',bilan1,bilan11
	1251	c print*,'Bilan2:',bilan2,bilan12
	1252	c print*,'Bilan3:',bilan3,bilan13
	1253	c print*,'Bilan5:',bilan5,bilan15
	1254
	1255	dice1=0.
	1256	dice2=0.
	1257	dice3=0.
	1258	dice4=0.
	1259	dice1=(bilan11-bilan1)rhoi_ch4 !glace 1 m^3.m^-2 kg.m^-3 pourquoi ????
	1260	dice2=(bilan12-bilan2)*rhoi_c2h6 !glace 2
	1261	dice3=(bilan13-bilan3)*rhoi_c2h2 !glace 3
	1262	dice4=(bilan14-bilan4)*rhol !noyaux
	1263
	1264	return
	1265	end
	1266
	1267
	1268	subroutine coagul2(ihor)
	1269
	1270	*********************************************************
	1271	* ce programme calcule la nouvelle concentration dans *
	1272	* le a ieme intervalle de rayon, a l'altitude h, a *
	1273	* l'instant t+dt *
	1274	*********************************************************
	1275
	1276	use dimphy
	1277	IMPLICIT NONE
	1278	#include "dimensions.h"
	1279	#include "microtab.h"
	1280	#include "varmuphy.h"
	1281
	1282	* declaration des blocs communs
	1283	*------------------------------
	1284
	1285	common/x2ctps/li,lf,dt
	1286	common/x2con/c,c1,c2,c3,caer
	1287	common/cldpart/rad,mmu
	1288
	1289	* declaration des variables
	1290	* --------------------------
	1291
	1292	integer li,lf
	1293	real dt
	1294	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
	1295	& c3(nz,nrad,2)
	1296	real caer(nz,nrad,2)
	1297	real rad(nz,nrad),mmu(nz,nrad)
	1298
	1299	* declaration des variables propres au ss-programme
	1300	* -------------------------------------------------
	1301
	1302	integer h,a,ihor,i
	1303	real pr,pe,eta1,eta2
	1304	real sum11,sum12,sum13,sum21,sum22,sum23
	1305	real rfb,rfx,rpr
	1306
	1307	* traitement
	1308	* ----------
	1309
	1310
	1311	sum11=0.
	1312	sum12=0.
	1313	sum13=0.
	1314	sum21=0.
	1315	sum22=0.
	1316	sum23=0.
	1317
	1318	do 10 h=nztop,nz
	1319	do 11 a=1,nrad ! boucle aerosol secs
	1320	call pertpro2(ihor,h,a,pe,pr)
	1321
	1322	if(1+dtpe .gt. 1.1) print,a,1+dt*pe,' scav'
	1323	caer(h,a,lf)=(caer(h,a,li)+prdt)/(1+dtpe)
	1324	c(h,a,lf)=c(h,a,li)
	1325	c1(h,a,lf)=c1(h,a,li)
	1326	c2(h,a,lf)=c2(h,a,li)
	1327	c
	1328	c eta1=0.
	1329	c eta2=0.
	1330	c if(c(h,a,li) .ne. 0.) eta1=c1(h,a,li)/c(h,a,li)
	1331	c if(c(h,a,li) .ne. 0.) eta2=c2(h,a,li)/c(h,a,li)
	1332	c c(h,a,lf) =( c(h,a,li)+prdt)/(1+dtpe)
	1333	c c1(h,a,lf)=(c1(h,a,li)+eta1prdt)/(1+dt*pe)
	1334	c c2(h,a,lf)=(c2(h,a,li)+eta2prdt)/(1+dt*pe)
	1335	c
	1336	c
	1337	c sum11=sum11+c(h,a,li)
	1338	c sum12=sum12+c1(h,a,li)
	1339	c sum13=sum13+c2(h,a,li)
	1340	c
	1341	c sum21=sum21+ c(h,a,lf)
	1342	c sum22=sum22+c1(h,a,lf)
	1343	c sum23=sum23+c2(h,a,lf)
	1344	c
	1345	c
	1346	11 continue
	1347	10 continue
	1348
	1349
	1350	if (nztop.ne.1) then
	1351	do 12 h=1,nztop-1
	1352	do 12 a=1,nrad
	1353	c(h,a,lf)=c(h,a,li)
	1354	c1(h,a,lf)=c1(h,a,li)
	1355	c2(h,a,lf)=c2(h,a,li)
	1356	caer(h,a,lf)=caer(h,a,li)
	1357	12 continue
	1358	endif
	1359
	1360	return
	1361	end
	1362
	1363
	1364	*__________________________________________________________________________
	1365
	1366	subroutine calcoag2
	1367
	1368	***************************************************************
	1369	* *
	1370	* Ce programme calcule les coefficients de collection d'une *
	1371	* particule de rayon x avec une particule de rayon b a une *
	1372	* altitude donnee h *
	1373	***************************************************************
	1374
	1375	* declaration des blocs communs
	1376	*------------------------------
	1377	use dimphy
	1378	IMPLICIT NONE
	1379	#include "dimensions.h"
	1380	#include "microtab.h"
	1381	#include "varmuphy.h"
	1382
	1383	common/x2ctps/li,lf,dt
	1384	common/x2con/c,c1,c2,c3,caer
	1385	common/cldpart/rad,mmu
	1386	common/x2coag/k
	1387	common/x2frac/rfg,dfg
	1388
	1389	* declaration des variables
	1390	* --------------------------
	1391
	1392	integer li,lf
	1393	real dt
	1394	real knu2,nud2,k(nz,nrad,nrad)
	1395	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
	1396	& c3(nz,nrad,2)
	1397	real caer(nz,nrad,2)
	1398	real rad(nz,nrad),mmu(nz,nrad)
	1399	real rfg(nz),dfg(nz,nrad)
	1400
	1401
	1402	* declaration des variables propres au ss-programme
	1403	* -------------------------------------------------
	1404
	1405	integer h,b,x,ihor,i
	1406	real nua,lambb,lambx,knb,knx,alphab,alphax,d,e,f,kcg
	1407	real db,dx,rm,dm,deltab,deltax,del,g,beta,gx,gb
	1408	real rfb,rfx,rpr
	1409	real*8 ne,qe,epso
	1410	real*8 corelec,yy
	1411
	1412	real kco,vx,vb,vitesse2,sto,ee,a,dd,bb,p0,t0,l0,ccol
	1413	real st(37),ef(37)
	1414	real vitesse,knu
	1415	external vitesse,knu
	1416
	1417
	1418	* initialisation
	1419	* --------------
	1420
	1421
	1422
	1423	* -nombres de STOCKES
	1424
	1425	data(st(i),i=1,37)/1.35,1.5,1.65,1.85,2.05,2.25,2.5,2.8,3.1,
	1426	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,
	1427	s 10.9,11.1,13.5,15.3,17.25,20.5,24.5,30.4,39.3,48,57,86.,
	1428	s 187.,600./
	1429
	1430	* -coef. d'efficacite de collection
	1431
	1432	ef(1)=3.75
	1433	ef(2)=8.75
	1434	do 11 i=3,37
	1435	ef(i)=ef(i-1)+2.5
	1436	11 continue
	1437
	1438	do 2 i=1,37
	1439	ef(i)=ef(i)*1e-2
	1440	2 continue
	1441
	1442	qe=1.6e-19
	1443	ne=-30e+6
	1444	epso=1e-9/(36*pi)
	1445
	1446	d=1.257
	1447	e=0.4
	1448	f=-1.1
	1449
	1450	* iteration sur z
	1451
	1452	do 1 h=1,nz
	1453
	1454	nua=nud2(h,1)
	1455
	1456	* iteration sur les rayons
	1457
	1458	do 1 b=1,nrad ! boucle aerosols secs : indice ''
	1459
	1460	knb=knu(h,b,1) ! knu et vitesse se trouvent dans brume.F
	1461	vb=vitesse(h,b,1) ! ils concernent les aerosols
	1462
	1463	do 1 x=1,nrad !boucles gouttes : indice '2'
	1464
	1465	knx=knu2(h,x,1)
	1466	vx=vitesse2(h,x,1)
	1467
	1468	COAGULATION **************************************************
	1469	--------------**************************************************
	1470	* calcul du terme correcteur 'slip-flow'
	1471
	1472	alphab=d+e*exp(f/knb)
	1473	alphax=d+e*exp(f/knx)
	1474
	1475	* calcul du coefficient de diffusion
	1476
	1477
	1478	rfb=(r_e(b)*(3./df(b))) ((rf(b)) **(1.-3./df(b)))
	1479	rfx=(rad(h,x)**
	1480	& (3./dfg(h,x)))((rfg(x))*(1.-3./dfg(h,x)))
	1481
	1482	db=kbzt(h)(1+alphabknb)/(6pinuarfb)
	1483	dx=kbzt(h)(1+alphaxknx)/(6pinuarfx)
	1484
	1485	* calcul du coefficient de coagulation
	1486
	1487	rpr=rfb+rfx
	1488	kcg=4pirpr*(db+dx)
	1489
	1490	* calcul de la vitesse thermique
	1491
	1492	gb=sqrt(6kbzt(h)/(rholpi2r_e(b)**3))
	1493	gx=sqrt(6kbzt(h)/(rholpi2rad(h,x)**3))
	1494
	1495	* calcul du libre parcours apparent des aerosols
	1496
	1497	lambb=8db/(pigb)
	1498	lambx=8dx/(pigx)
	1499
	1500	*calcul du terme correcteur beta
	1501
	1502	rm=rpr/2.
	1503	dm=(dx+db)/2.
	1504	g=sqrt(gx2+gb2)
	1505	deltab=(((2rfb+lambb)3-(4rfb2+lambb2)**1.5)
	1506	s /(6rfblambb)-2rfb)sqrt(2.)
	1507	deltax=(((2rfx+lambx)3-(4rfx2+lambx2)**1.5)
	1508	s /(6rfxlambx)-2rfx)sqrt(2.)
	1509	del=sqrt(deltab2+deltax2)
	1510	beta=1/((rm/(rm+del/2))+(4dm/(grm)))
	1511
	1512	* calcul du coefficient de coagulation corrige
	1513
	1514	kcg=kcg*beta
	1515
	1516	COALESCENCE ************************************************
	1517	-------------************************************************
	1518
	1519	kco=0.
	1520
	1521	if ( b.eq. x) continue ! goto 9
	1522
	1523
	1524	* calcul du nombre de Stockes de la petite particule
	1525
	1526	sto=2rholrfx*2abs(vx-vb)/(9nuarfb)
	1527
	1528	* calcul du coef. de Cunningham-Millikan
	1529
	1530	a=1.246
	1531	bb=0.42
	1532	dd=0.87
	1533	l0=0.653e-7
	1534	p0=101325.
	1535	t0=288.
	1536
	1537	ee=1+
	1538	& (l0t(h)p0(a+bbexp(-ddrfxt0p(h)/(l0t(h)*p0))))
	1539	s /(rfxt0p(h))
	1540
	1541	* calcul du nombre de Stockes corrige
	1542
	1543	sto=sto*ee
	1544
	1545	if (sto .le. 1.2) goto 9
	1546
	1547	if (sto .ge. 600.) then
	1548	ccol=1.
	1549	goto 8
	1550	endif
	1551
	1552	* recherche du coefficient de collection
	1553
	1554	do 3 i=1,37
	1555	if (sto .gt. st(i)) then
	1556	goto 3
	1557	endif
	1558	if (sto .eq. st(i)) then
	1559	ccol=ef(i+1)
	1560	else
	1561	ccol=ef(i)
	1562	endif
	1563	goto 8
	1564	3 continue
	1565
	1566	* calcul du coefficient de coalescence
	1567
	1568	8 kco=pi(rfb+rfx)2ccol*abs(vb-vx)
	1569
	1570	9 continue
	1571
	1572	CORRECTION ELECTRICITE *****************************
	1573	------------------------****************************
	1574
	1575	c yy=1.d0ne2r(x)r(b)qe**2
	1576	c & /(1.d0kbzt(h)(r(b)+r(x))4piepso)
	1577	c corelec=0.
	1578	c if (yy.lt.50.) corelec=yy/(exp(yy)-1.)
	1579
	1580	corelec=1.
	1581
	1582	c b=aerosol
	1583	c x=gouttes
	1584
	1585	k(h,b,x)=(kcg+kco)*corelec
	1586
	1587	c
	1588	c ATTENTION, IL N'Y A PLUS DE SYMETRIE...
	1589	c k(ihor,h,x,b)=k(ihor,h,b,x)
	1590	c
	1591	c
	1592
	1593
	1594	1 continue
	1595	return
	1596	end
	1597
	1598	*______________________________________________________________________
	1599
	1600	subroutine pertpro2(ihor,h,a,l_,pr_)
	1601
	1602	*****************************************************************************
	1603	* *
	1604	* ce programme permet le calcul du terme de production (pr) et de perte (l)*
	1605	* pour le phenomene de coagulation *
	1606	* dans le a ieme intervalle de rayon a une altitude h *
	1607	****************************************************************************
	1608
	1609	* declaration des blocs communs
	1610	*------------------------------
	1611	use dimphy
	1612	IMPLICIT NONE
	1613	#include "dimensions.h"
	1614	#include "microtab.h"
	1615	#include "varmuphy.h"
	1616
	1617	common/x2ctps/li,lf,dt
	1618	common/x2con/c,c1,c2,c3,caer
	1619	common/cldpart/rad,mmu
	1620	common/x2coag/k
	1621
	1622
	1623	* declaration des variables
	1624	* --------------------------
	1625
	1626	integer li,lf
	1627	real dt
	1628	real dr(nrad),dv(nrad)
	1629	real k(nz,nrad,nrad)
	1630	real c(nz,nrad,2), c1(nz,nrad,2), c2(nz,nrad,2),
	1631	& c3(nz,nrad,2)
	1632	real caer(nz,nrad,2)
	1633	real rad(nz,nrad),mmu(nz,nrad)
	1634
	1635
	1636	* declaration des variables propres au ss-programme
	1637	* -------------------------------------------------
	1638
	1639	integer h,b,a,x,ihor,i
	1640	real*8 pr,ss,s,l
	1641	real pr_,l_,vol,del
	1642
	1643	* traitement
	1644	* -----------
	1645
	1646
	1647
	1648	* production
	1649	*+++++++++++++
	1650	s=0.d0
	1651	ss=0.d0
	1652	pr=0.d0
	1653
	1654	c Pas de production d'aerosols par scavenging !!!
	1655
	1656	pr=0.d0
	1657
	1658	* perte
	1659	*- - - - -
	1660
	1661	l=0.d0
	1662
	1663
	1664	do 10 x=1,nrad ! boucle sur les gouttes
	1665	l=l+k(h,a,x)*c(h,x,li)
	1666	if(l.ne.0.d0) then
	1667	print*,a,x,k(h,a,x),c(h,x,li),l,' : detail coal'
	1668	endif
	1669	10 continue
	1670
	1671
	1672	#ifdef CRAY
	1673	l_=l
	1674	pr_=pr
	1675	#else
	1676	l_=sngl(l)
	1677	pr_=sngl(pr)
	1678	#endif
	1679
	1680	return
	1681
	1682	end
	1683

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

Download in other formats:

Original Format