Context Navigation

source: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c @ 1242

Last change on this file since 1242 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 32.3 KB

Rev	Line
[3]	1	/* gptitan: photochimie */
	2	/* GCCM */
	3
	4	/* tout est passe en simple precision */
	5	/* sauf pour l'inversion de la matrice */
	6
	7	/* nitriles et hydrocarbures separes pour l'inversion */
	8
	9	/* flux variable au sommet */
	10
	11	#include "titan.h"
	12
[1057]	13	void gptitan_(
[3]	14	double RA, double TEMP, double *NB,
[1126]	15	char CORPS[][10], double Y[][NLEV],
	16	double FIN, int LAT, double *MASS, double MD[][NLEV],
	17	double KEDD, double botCH4, double KRATE[][NLEV],
[3]	18	int reactif[][5], int nom_prod, int nom_perte,
	19	int prod[][200], int perte[][200][2], int aerprod, int utilaer,
	20	double MAER[][NLEV], double PRODAER[][NLEV],
	21	double CSN[][NLEV], double CSH[][NLEV],
	22	int htoh2, double surfhaze)
	23	{
	24	char outlog[100],corps[100][10];
[1126]	25	int i,j,k,l;
[3]	26	int ireac,ncom1,ncom2;
[1126]	27	double *a,b,*c;
	28	double fl,fp,mu,jac,ym1,*f;
	29	double fluxCH4;
	30	double conv,delta,deltamax;
	31	double cm,cp,dim,dip,dm,dp,dym,dyp,km,kp,r,dra,dram,drap;
	32	double np,nm,s,test,time,ts,v,dv;
[3]	33	char str2[15];
	34	FILE *out;
	35
	36	/* va avec htoh2 */
	37	double dyh,dyh2;
	38
	39	/* va avec aer */
	40	double dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
	41	double k_dep,faer;
	42	double productaer,csurn,csurh,mmolaer;
	43
	44	if( (*aerprod) == 1 )
	45	{
	46	k_dep = dm2d( 1, 5, 1, 3 ); /* k en s-1, reactions d'initiation */
	47	faer = dm2d( 1, 5, 1, 3 ); /* fraction de chaque compose */
	48	productaer = dm1d( 0, 3 ); /* local production rate by pathways */
	49	mmolaer = dm1d( 0, 3 ); /* local molar mass by pathways */
	50	csurn = dm1d( 0, 3 ); /* local C/N by pathways */
	51	csurh = dm1d( 0, 3 ); /* local C/H by pathways */
	52	}
	53
	54	/* DEBUG */
[1126]	55	printf("CHIMIE: lat=%d\n",(*LAT)+1);
[3]	56	/**/
	57
	58	for( i = 0; i <= NC; i++)
	59	{
	60	strcpy( corps[i], CORPS[i] );
	61	corps[i][strcspn(CORPS[i], " ")] = '\0';
	62	}
	63
	64	strcpy( outlog, "chimietitan" );
	65	strcat( outlog, ".log" );
[1126]	66	out = fopen( outlog, "w" );
	67	fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
	68	fclose( out );
	69
[3]	70	deltamax = 1.e5;
[1126]	71	test = 1.0e-15;
[3]	72
[1126]	73	/* valeur de r:
	74	r = g0 R0^2 / R * 2 * 1E-3
	75	avec g0 en cm/s2, R0 en km, mu et mass en g
	76	*/
	77	r = 21.595656e0;
	78
[3]	79	/* DEBUG
	80	out = fopen( outlog, "a" );
[1126]	81	fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
[3]	82	fclose( out );
	83	*/
	84	fl = dm1d( 0, NC-1 );
	85	fp = dm1d( 0, NC-1 );
	86	mu = dm1d( 0, NLEV-1 );
[1126]	87	ym1 = dm2d( 0, NC-1, 0, NLEV-1 );
	88	f = dm2d( 0, NC-1, 0, NLEV-1 );
[3]	89	jac = dm2d( 0, NC-1, 0, NC-1 );
[1126]	90	c = dm2d( 0, NLEV-1, 0, NC-1 );
	91	a = dm3d( 0, NLEV-1, 0, NC-1, 0, NC-1 );
	92	b = dm3d( 0, NLEV-1, 0, NC-1, 1, 2 );
[3]	93
	94	/* DEBUG */
	95	/*
	96	out = fopen( "err.log", "a" );
	97	fprintf( out,"%s\n", );
	98	fclose( out );
	99	*/
	100
	101	/* initialisation mu, CH4 au sol */
	102
	103	for( j = 0; j <= NLEV-1; j++ )
	104	{
	105	mu[j] = 0.0e0;
	106	for( i = 0; i <= ST-1; i++ )
	107	{
	108	if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
	109	{
	110	fluxCH4 = (*botCH4 - Y[i][j]);
	111	Y[i][j] = *botCH4;
	112	}
	113	mu[j] += ( MASS[i] * Y[i][j] );
	114	}
	115	}
	116
[1126]	117	/* initialisation compo avant calcul */
	118	for( j = NLEV-1; j >= 0; j-- )
	119	for( i = 0; i <= ST-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
	120
	121	/*
	122	==========================================================================
	123	STRATEGIE:
	124	INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
	125	PUIS INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
	126	==========================================================================
	127
	128	PREMIERE ETAPE:
	129	===============
	130	INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
	131	===============
	132	*/
	133
[3]	134	/* ****************** */
[1126]	135	/* Time loop: */
[3]	136	/* ****************** */
	137
[1126]	138	time = ts = 0.0e0;
	139	delta = 1.e-3;
[3]	140
[1126]	141	while( time < (*FIN) )
	142	{
	143
	144
[3]	145	/* DEBUG
[1126]	146	for( j = NLEV-1; j >= NLD; j-- )
	147	{
[3]	148	out = fopen( outlog, "a" );
	149	fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
	150	fclose( out );
	151
	152	out = fopen( "profils.log", "a" );
	153	fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
	154	for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
	155	for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
	156	fclose( out );
[1126]	157	}
	158	exit(0);
	159	*/
[3]	160
	161
[1126]	162	/* ------------------------------ */
	163	/* Calculs variations et jacobien */
	164	/* ------------------------------ */
	165
	166	for( j = NLEV-1; j >= NLD; j-- )
	167	{
	168
	169	/* init of step */
	170	/* ------------ */
	171	for( i = 0; i <= ST-1; i++ )
	172	{
	173	fp[i] = fl[i] = 0.0e0;
	174	for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
	175	}
	176
	177	/* Chimie */
	178	/* ------ */
	179
	180	/* productions et pertes chimiques */
	181	for( i = 0; i <= ST-1; i++ )
	182	{
	183	Y[i][j] = max(Y[i][j],1.e-30); /* minimum */
	184
	185	for( l = 0; l <= nom_prod[i]-1; l++ ) /* Production term */
	186	{
	187	ireac = prod[i][l]; /* Number of the reaction involves. */
	188	ncom1 = reactif[ireac][0]; /* First compound which reacts. */
	189	if( reactif[ireac][1] == NC ) /* Photodissociation or relaxation */
	190	{
	191	jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
	192	fp[i] += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
	193	}
	194	else /* General case. */
	195	{
	196	ncom2 = reactif[ireac][1]; /* Second compound which reacts. */
	197	jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] ); /* Jacobian compound #1. */
	198	jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] ); /* Jacobian compound #2. */
	199	fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
	200	}
	201	}
	202
	203	for( l = 0; l <= nom_perte[i]-1; l++ ) /* Loss term. */
	204	{
	205	ireac = perte[i][l][0]; /* Reaction number. */
	206	ncom2 = perte[i][l][1]; /* Compound #2 reacts. */
	207	if( reactif[ireac][1] == NC ) /* Photodissociation or relaxation */
	208	{
	209	jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
	210	fl[i] += ( KRATE[ireac][j] * NB[j] );
	211	}
	212	else /* General case. */
	213	{
	214	jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] ); /* Jacobian compound #1. */
	215	jac[i][i] -= ( KRATE[ireac][j] * Y[ncom2][j] ); /* Jacobien compound #2. */
	216	fl[i] += ( KRATE[ireac][j] * Y[ncom2][j] ); /* Loss term. */
	217	}
	218	}
	219	}
	220
	221
	222	/* Aerosols */
	223	/* -------- */
	224	if( (*aerprod) == 1 )
	225	{
	226	aer(corps,TEMP,NB,Y,&j,k_dep,faer,
	227	&dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
	228	mmolaer,productaer,csurn,csurh);
	229
	230	for( i = 0; i <= 3; i++ )
	231	{
	232	PRODAER[i][j] = productaer[i];
	233	MAER[i][j] = mmolaer[i];
	234	CSN[i][j] = csurn[i];
	235	CSH[i][j] = csurh[i];
	236	}
	237	/* DEBUG
	238	printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
	239	if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
	240	printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
	241	fp[utilaer[2]],dyc2h2*NB[j]);
	242	if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
	243	printf("fp(%s) =%e; dyhcn =%e\n",corps[utilaer[5]],
	244	fp[utilaer[5]],dyhcn*NB[j]);
	245	if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
	246	printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
	247	fp[utilaer[6]],dyhc3n*NB[j]);
	248	if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
	249	printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
	250	fp[utilaer[13]],dynccn*NB[j]);
	251	if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
	252	printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
	253	fp[utilaer[14]],dych3cn*NB[j]);
	254	if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
	255	printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
	256	fp[utilaer[15]],dyc2h3cn*NB[j]);
[3]	257	*/
	258
[1126]	259	fp[utilaer[2]] -= ( dyc2h2 * NB[j] );
	260	fp[utilaer[5]] -= ( dyhcn * NB[j] );
	261	fp[utilaer[6]] -= ( dyhc3n * NB[j] );
	262	fp[utilaer[13]]-= ( dynccn * NB[j] );
	263	fp[utilaer[14]]-= ( dych3cn * NB[j] );
	264	fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
	265	if( Y[utilaer[2]][j] != 0.0 )
	266	jac[utilaer[2]][utilaer[2]] -= ( dyc2h2 * NB[j] / Y[utilaer[2]][j] );
	267	if( Y[utilaer[5]][j] != 0.0 )
	268	jac[utilaer[5]][utilaer[5]] -= ( dyhcn * NB[j] / Y[utilaer[5]][j] );
	269	if( Y[utilaer[6]][j] != 0.0 )
	270	jac[utilaer[6]][utilaer[6]] -= ( dyhc3n * NB[j] / Y[utilaer[6]][j] );
	271	if( Y[utilaer[13]][j] != 0.0 )
	272	jac[utilaer[13]][utilaer[13]] -= ( dynccn * NB[j] / Y[utilaer[13]][j] );
	273	if( Y[utilaer[14]][j] != 0.0 )
	274	jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
	275	if( Y[utilaer[15]][j] != 0.0 )
	276	jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
	277	}
	278
	279
	280	/* H -> H2 on haze particles */
	281	/* ------------------------- */
	282	if( (*htoh2) == 1 )
	283	{
	284	heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
	285	/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
	286
	287	/* DEBUG
	288	printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
	289	if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
	290	printf("fp(%s) = %e; dyh = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
	291	if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
	292	printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
	293	*/
	294
	295	fp[utilaer[0]] += ( dyh * NB[j] );
	296	/* pourquoi pas 2 ?? cf gptit dans 2da... /
	297
	298	fp[utilaer[1]] += ( dyh2 * NB[j] );
	299	if( Y[utilaer[0]][j] != 0.0 )
	300	jac[utilaer[0]][utilaer[0]] += ( dyh * NB[j] / Y[utilaer[0]][j] );
	301	/* pourquoi pas 2 ?? cf gptit dans 2da... /
	302	}
	303
	304
	305	/* Backup jacobian level j. */
	306	/* ------------------------ */
	307	for( i = 0; i <= ST-1; i++ )
	308	for( k = 0; k <= ST-1; k++ )
	309	a[j][i][k] = jac[i][k];
	310
	311
	312	/* Diffusion verticale et flux exterieurs */
	313	/* -------------------------------------- */
	314
	315	/*
	316	pour dy/dr, dr doit etre en cm...
	317	pareil pour dphi/dr
	318	*/
	319	for( i = 0; i <= ST-1; i++ )
	320	{
	321
	322	/* First level. */
	323	if( j == NLD )
	324	{
	325	v = dv = 0.0e0;
	326	dra = RA[j+1]-RA[j];
	327
	328	cp = (NB[j+1]+NB[j])/2.; /* Mean total concentration. */
	329	dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
	330	pow( RA[j+1], 2.0e0 ); /* Delta i,j level +1. */
	331	dp = (MD[i][j]+MD[i][j+1])/2.; /* Mean molecular diffusion. */
	332	dyp = (Y[i][j+1]-Y[i][j])/(RA[j+2]-RA[j])2.e-5; / Delta y level +1. */
	333	kp = (KEDD[j+1]+KEDD[j])/2.; /* Mean eddy diffusion. */
	334	/* div phi. */
	335	f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
	336	+ kp * dyp )
	337	* (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.))
	338	+ fp[i] - Y[i][j]*fl[i] + v;
	339	/* dphi / dy this level. */
	340	a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
	341	- 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
	342	* (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) + dv );
	343	/* dphi / dy level +1. */
	344	c[j][i] = -THETA * delta
	345	* cp * ( dp * 0.5e0 * dip
	346	+ 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
	347	* (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.));
	348	}
	349	/* Last level. */
	350	else if( j == NLEV-1 )
	351	{
	352	v = dv = 0.0e0;
	353	dra = RA[NLEV-1]-RA[NLEV-2];
	354
	355	/* Jeans escape */
	356	if( strcmp(corps[i], "H") == 0 )
	357	{
	358	dv = top_H * NB[NLEV-1]
	359	* (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
	360	v = dv * Y[i][NLEV-1];
	361	}
	362	if( strcmp(corps[i], "H2") == 0 )
	363	{
	364	dv = top_H2 * NB[NLEV-1]
	365	* (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
	366	v = dv * Y[i][NLEV-1];
	367	}
	368	/* Input flux for N(4S) */
	369	if( strcmp(corps[i], "N4S") == 0 )
	370	v = top_N4S
	371	* (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
	372
	373	cm = (NB[NLEV-1]+NB[NLEV-2])/2.; /* Mean total concentration. */
	374	dim = r * (MASS[i]-(mu[NLEV-1]+mu[NLEV-2])/2.)
	375	/ (TEMP[NLEV-1]+TEMP[NLEV-2])
	376	/ pow( RA[NLEV-1], 2.0e0 ); /* Delta i,j level -1. */
	377	dm = (MD[i][NLEV-1]+MD[i][NLEV-2])/2.; /* Mean molecular diffusion. */
	378	dym = (Y[i][NLEV-1]-Y[i][NLEV-2])/dra1.e-5; / Delta y level -1. */
	379	km = (KEDD[NLEV-1]+KEDD[NLEV-2])/2.; /* Mean eddy diffusion. */
	380	/* div phi. */
	381	f[i][NLEV-1] = fp[i] - Y[i][NLEV-1]*fl[i] - v
	382	- cm * ( dm * ( (Y[i][NLEV-1]+Y[i][NLEV-2])/2. * dim + dym )
	383	+ km * dym )
	384	* (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
	385	/* dphi / dy this level */
	386	a[NLEV-1][i][i] -= ( cm * ( dm * 0.5e0 * dim
	387	+ 1.e-5/dra * (dm + km ) )
	388	* (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)) + dv );
	389	/* dphi / dy level -1. */
	390	b[NLEV-1][i][2] = THETA * delta
	391	* cm * ( dm * 0.5e0 * dim
	392	- 1.e-5/dra * (dm + km ) )
	393	* (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
	394	}
	395	else
	396	{
	397	v = dv = 0.0e0;
	398	dram=(RA[j+1]-RA[j-1])/2.;
	399	if (j<NLEV-2)
	400	drap=(RA[j+1]-RA[j-1])/2.;
	401	else
	402	drap=dram;
	403
	404	cm = (NB[j]+NB[j-1])/2.; /* Mean concentration level -1. */
	405	cp = (NB[j]+NB[j+1])/2.; /* Mean concentration level +1. */
	406	dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
	407	pow( RA[j+1], 2.0e0 ); /* Delta i,j level +1. */
	408	dim = r * (MASS[i]-(mu[j]+mu[j-1])/2.) / (TEMP[j]+TEMP[j-1]) /
	409	pow( RA[j], 2.0e0 ); /* Delta i,j level -1. */
	410	dm = (MD[i][j-1]+MD[i][j])/2.; /* Mean molecular diffusion level -1. */
	411	dp = (MD[i][j+1]+MD[i][j])/2.; /* Mean molecular diffusion level +1. */
	412	dym = (Y[i][j]-Y[i][j-1])/dram1.e-5; / Delta y level -1. */
	413	dyp = (Y[i][j+1]-Y[i][j])/drap1.e-5; / Delta y level +1. */
	414	km = (KEDD[j]+KEDD[j-1])/2.; /* Mean eddy diffusion level -1. */
	415	kp = (KEDD[j]+KEDD[j+1])/2.; /* Mean eddy diffusion level +1. */
	416	/* div phi. */
	417	f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
	418	+ kp * dyp )
	419	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
	420	- cm * ( dm * ( (Y[i][j]+Y[i][j-1])/2. * dim + dym )
	421	+ km * dym )
	422	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.))
	423	+ fp[i] - fl[i] * Y[i][j] + v;
	424	/* dphi / dy this level */
	425	a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
	426	- 1.e-5/drap * (dp + kp) )
	427	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
	428	- cm * ( dm * 0.5e0 * dim
	429	+ 1.e-5/dram * (dm + km ) )
	430	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) );
	431	/* dphi / dy level -1. */
	432	b[j][i][2] = THETA * delta
	433	* cm * ( dm * 0.5e0 * dim
	434	- 1.e-5/dram * (dm + km ) )
	435	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.));
	436	/* dphi / dy level +1. */
	437	c[j][i] = -THETA * delta
	438	* cp * ( dp * 0.5e0 * dip
	439	+ 1.e-5/drap * (dp + kp) )
	440	* (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.));
	441	}
	442	}
	443
	444
	445
	446	/* finition pour inversion */
	447	/* ----------------------- */
	448
	449	for( i = 0; i <= ST-1; i++ )
	450	{
	451	for( k = 0; k <= ST-1; k++ )
	452	{
	453	a[j][i][k] = ( -THETA delta ); /* Correction time step. */
	454	if( k == i ) a[j][k][k] += NB[j]; /* Correction diagonal. */
	455	}
	456	f[i][j] *= delta;
	457	}
	458
	459	}
	460
	461
	462	/* -------------------------------- */
	463	/* Inversion of matrix cf method LU */
	464	/* -------------------------------- */
	465
	466	for( j = NLD+1; j <= NLEV-1; j++ )
	467	{
	468	solve( a, j-1, 0, ST-1 );
	469	for( i = 0; i <= ST-1; i++ )
	470	{
	471	s = 0.0e0;
	472	for( k = 0; k <= ST-1; k++ )
	473	{
	474	a[j][i][k] -= ( b[j][i][2] * c[j-1][k] * a[j-1][i][k] );
	475	s += ( b[j][i][2] * f[k][j-1] * a[j-1][i][k] );
	476	}
	477	f[i][j] -= s;
	478	}
	479	}
	480	solve( a, NLEV-1, 0, ST-1 );
	481	for( j = NLEV-1; j >= NLD; j-- )
	482	{
	483	if( j != NLEV-1 )
	484	for( i = 0; i <= ST-1; i++ ) f[i][j] -= ( c[j][i] * b[j+1][i][1] );
	485	for( i = 0; i <= ST-1; i++ )
	486	{
	487	s = 0.0e0;
	488	for( k = 0; k <= ST-1; k++ ) s += ( a[j][i][k] * f[k][j] );
	489	b[j][i][1] = s;
	490	Y[i][j] += s;
	491	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
	492	}
	493	}
	494
	495	/* ------------------ */
	496	/* Tests et evolution */
	497	/* ------------------ */
	498
	499	/* Calcul deviation */
	500	/* ---------------- */
	501
	502	for( j = NLD; j <= NLEV-1; j++ )
	503	for( i = 0; i <= ST-1; i++ )
	504	if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
	505	{
	506	conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
	507	if( conv > ts )
	508	{
	509	/*
	510	if( conv >= 0.1 )
	511	{
	512	out = fopen( outlog, "a" );
	513	fprintf( out, "Lat no %d;", (*LAT)+1);
	514	fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
	515	fclose( out );
	516	}
	517	*/
	518	ts = conv;
	519	}
	520	}
	521
	522	/* test deviation */
	523	/* -------------- */
	524
	525	if( ts < 0.1e0 )
	526	{
	527	for( i = 0; i <= ST-1; i++ )
	528	for( j = NLD; j <= NLEV-1; j++ )
	529	if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
	530	{
	531	out = fopen( outlog, "a" );
	532	fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
	533	corps[i], ym1[i], Y[i][j], j );
	534	for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
	535	fclose( out );
	536	exit(0);
	537	// Y[i][j] = 1.e-20;
	538	}
	539	for( j = NLD; j <= NLEV-1; j++ )
	540	for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
	541	time += delta;
	542	if( ts < 1.00e-5 ) delta *= 1.0e2;
	543	if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
	544	if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
	545	if( ( ts > 1.00e-3 ) && ( ts < 5.0e-3 ) ) delta *= 3.0e0;
	546	if( ( ts > 5.00e-3 ) && ( ts < 0.01e0 ) ) delta *= 1.5e0;
	547	if( ( ts > 0.010e0 ) && ( ts < 0.03e0 ) ) delta *= 1.2e0;
	548	if( ( ts > 0.030e0 ) && ( ts < 0.05e0 ) ) delta *= 1.1e0;
	549
	550	// if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
	551	// if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
	552
	553	delta = min( deltamax, delta );
	554	}
	555	else
	556	{
	557	for( j = NLD; j <= NLEV-1; j++ )
	558	for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
	559
	560	if( ts > 0.8 ) delta *= 1.e-6;
	561	if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
	562	if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
	563	if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
	564	if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
	565	if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
	566	}
	567	ts = 0.0e0;
	568
	569	out = fopen( outlog, "a" );
	570	fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
	571	fclose( out );
	572
	573	}
	574	/* **************** */
	575	/* end of time loop */
	576	/* **************** */
	577
	578	/*
	579	==========================================================================
	580
	581	SECONDE ETAPE:
	582	===============
	583	INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
	584	===============
	585	*/
	586	if( NLD != 0 )
	587	for( j = NLD-1; j >= 0; j-- )
	588	{
[3]	589	time = ts = 0.0e0;
	590	delta = 1.e-3;
	591
	592	/* ++++++++++++ */
	593	/* time loop. */
	594	/* ++++++++++++ */
	595
	596	while( time < (*FIN) )
	597	{
	598
[1126]	599	/* init of step */
	600	/* ------------ */
	601	for( i = 0; i <= ST-1; i++ )
	602	{
	603	fp[i] = fl[i] = 0.0e0;
	604	for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
	605	}
[3]	606
[1126]	607	/* Chimie */
	608	/* ------ */
	609
	610	/* productions et pertes chimiques */
	611	for( i = 0; i <= ST-1; i++ )
[3]	612	{
	613	Y[i][j] = max(Y[i][j],1.e-30); /* minimum */
	614
	615	for( l = 0; l <= nom_prod[i]-1; l++ ) /* Production term */
	616	{
	617	ireac = prod[i][l]; /* Number of the reaction involves. */
	618	ncom1 = reactif[ireac][0]; /* First compound which reacts. */
	619	if( reactif[ireac][1] == NC ) /* Photodissociation or relaxation */
	620	{
	621	jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
	622	fp[i] += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
	623	}
	624	else /* General case. */
	625	{
	626	ncom2 = reactif[ireac][1]; /* Second compound which reacts. */
	627	jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] ); /* Jacobian compound #1. */
	628	jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] ); /* Jacobian compound #2. */
	629	fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
	630	}
	631	}
	632
	633	for( l = 0; l <= nom_perte[i]-1; l++ ) /* Loss term. */
	634	{
	635	ireac = perte[i][l][0]; /* Reaction number. */
	636	ncom2 = perte[i][l][1]; /* Compound #2 reacts. */
	637	if( reactif[ireac][1] == NC ) /* Photodissociation or relaxation */
	638	{
	639	jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
	640	fl[i] += ( KRATE[ireac][j] * NB[j] );
	641	}
	642	else /* General case. */
	643	{
	644	jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] ); /* Jacobian compound #1. */
	645	jac[i][i] -= ( KRATE[ireac][j] * Y[ncom2][j] ); /* Jacobien compound #2. */
	646	fl[i] += ( KRATE[ireac][j] * Y[ncom2][j] ); /* Loss term. */
	647	}
	648	}
	649	}
	650
[1126]	651
[3]	652	/* Aerosols */
	653	/* -------- */
	654	if( (*aerprod) == 1 )
	655	{
	656	aer(corps,TEMP,NB,Y,&j,k_dep,faer,
	657	&dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
	658	mmolaer,productaer,csurn,csurh);
	659
	660	for( i = 0; i <= 3; i++ )
	661	{
	662	PRODAER[i][j] = productaer[i];
	663	MAER[i][j] = mmolaer[i];
	664	CSN[i][j] = csurn[i];
	665	CSH[i][j] = csurh[i];
	666	}
	667	/* DEBUG
	668	printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
	669	if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
	670	printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
	671	fp[utilaer[2]],dyc2h2*NB[j]);
	672	if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
	673	printf("fp(%s) =%e; dyhcn =%e\n",corps[utilaer[5]],
	674	fp[utilaer[5]],dyhcn*NB[j]);
	675	if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
	676	printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
	677	fp[utilaer[6]],dyhc3n*NB[j]);
	678	if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
	679	printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
	680	fp[utilaer[13]],dynccn*NB[j]);
	681	if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
	682	printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
	683	fp[utilaer[14]],dych3cn*NB[j]);
	684	if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
	685	printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
	686	fp[utilaer[15]],dyc2h3cn*NB[j]);
	687	*/
	688
	689	fp[utilaer[2]] -= ( dyc2h2 * NB[j] );
	690	fp[utilaer[5]] -= ( dyhcn * NB[j] );
	691	fp[utilaer[6]] -= ( dyhc3n * NB[j] );
	692	fp[utilaer[13]]-= ( dynccn * NB[j] );
	693	fp[utilaer[14]]-= ( dych3cn * NB[j] );
	694	fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
	695	if( Y[utilaer[2]][j] != 0.0 )
	696	jac[utilaer[2]][utilaer[2]] -= ( dyc2h2 * NB[j] / Y[utilaer[2]][j] );
	697	if( Y[utilaer[5]][j] != 0.0 )
	698	jac[utilaer[5]][utilaer[5]] -= ( dyhcn * NB[j] / Y[utilaer[5]][j] );
	699	if( Y[utilaer[6]][j] != 0.0 )
	700	jac[utilaer[6]][utilaer[6]] -= ( dyhc3n * NB[j] / Y[utilaer[6]][j] );
	701	if( Y[utilaer[13]][j] != 0.0 )
	702	jac[utilaer[13]][utilaer[13]] -= ( dynccn * NB[j] / Y[utilaer[13]][j] );
	703	if( Y[utilaer[14]][j] != 0.0 )
	704	jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
	705	if( Y[utilaer[15]][j] != 0.0 )
	706	jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
	707	}
[1126]	708
	709
[3]	710	/* H -> H2 on haze particles */
	711	/* ------------------------- */
	712	if( (*htoh2) == 1 )
	713	{
	714	heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
	715	/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
	716
	717	/* DEBUG
	718	printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
	719	if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
	720	printf("fp(%s) = %e; dyh = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
	721	if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
	722	printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
	723	*/
	724
[1126]	725	fp[utilaer[0]] += ( dyh * NB[j] );
	726	/* pourquoi pas 2 ?? cf gptit dans 2da... /
	727
[3]	728	fp[utilaer[1]] += ( dyh2 * NB[j] );
	729	if( Y[utilaer[0]][j] != 0.0 )
	730	jac[utilaer[0]][utilaer[0]] += ( dyh * NB[j] / Y[utilaer[0]][j] );
[1126]	731	/* pourquoi pas 2 ?? cf gptit dans 2da... /
[3]	732	}
	733
[1126]	734
	735	/* Backup jacobian level j. */
	736	/* ------------------------ */
	737	for( i = 0; i <= ST-1; i++ )
[3]	738	{
	739	for( k = 0; k <= ST-1; k++ )
[1126]	740	a[j][i][k] = jac[i][k];
	741	f[i][j] = fp[i] - fl[i] * Y[i][j];
[3]	742	}
	743
[1126]	744
	745	/* finition pour inversion */
	746	/* ----------------------- */
	747
	748	for( i = 0; i <= ST-1; i++ )
[3]	749	{
[1126]	750	for( k = 0; k <= ST-1; k++ )
	751	{
	752	a[j][i][k] = ( -THETA delta ); /* Correction time step. */
	753	if( k == i ) a[j][k][k] += NB[j]; /* Correction diagonal. */
	754	}
	755	f[i][j] *= delta;
[3]	756	}
	757
[1126]	758
[3]	759	/* Inversion of matrix cf method LU */
	760	/* -------------------------------- */
	761
	762	/* inversion by blocs: */
	763	/* Hydrocarbons */
	764
[1126]	765	solve_b( a, f, j, 0, NHC-1 );
[3]	766	for( i = 0; i <= NHC-1; i++ )
	767	{
[1126]	768	Y[i][j] += f[i][j];
[3]	769	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
	770	}
	771
	772	/* Nitriles */
	773
[1126]	774	solve_b( a, f, j, NHC, ST-1 );
[3]	775	for( i = NHC+1; i <= ST-1; i++ )
	776	{
[1126]	777	Y[i][j] += f[i][j];
[3]	778	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
	779	}
	780
	781	/* end inversion by blocs: */
	782
	783	/* CH4 au sol */
	784	/* ---------- */
[1126]	785	for( i = 0; i <= ST-1; i++ )
	786	if( ( strcmp(corps[i], "CH4") == 0 ) && (j==0) && ( Y[i][0] < *botCH4 ) )
	787	{
	788	fluxCH4 += (*botCH4 - Y[i][0]);
	789	Y[i][0] = *botCH4;
	790	}
[3]	791
[1126]	792	/* ------------------ */
	793	/* Tests et evolution */
	794	/* ------------------ */
[3]	795
[1126]	796	/* Calcul deviation */
	797	/* ---------------- */
[3]	798
	799	for( i = 0; i <= ST-1; i++ )
	800	{
	801	test = 1.0e-15;
[1126]	802	if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
[3]	803	{
[1126]	804	conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
[3]	805
	806	if( conv > ts )
	807	{
	808	/*
	809	if( conv >= 0.1 )
	810	{
	811	out = fopen( outlog, "a" );
	812	fprintf( out, "Lat no %d; declin:%e;", (LAT)+1, (DECLIN) );
	813	fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
	814	fclose( out );
	815	}
	816	*/
	817	ts = conv;
	818	}
	819	}
	820	}
	821
[1126]	822	/* test deviation */
	823	/* -------------- */
	824
[3]	825	if( ts < 0.1e0 )
	826	{
	827	for( i = 0; i <= ST-1; i++ )
	828	if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
	829	{
	830	out = fopen( outlog, "a" );
	831	fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
[1126]	832	corps[i], ym1[i][j], Y[i][j], j );
	833	for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i][j],Y[i][k] );
[3]	834	fclose( out );
	835	// exit(0);
	836	Y[i][j] = 1.e-20;
	837	}
[1126]	838	for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
[3]	839	time += delta;
	840	if( ts < 1.00e-5 ) delta *= 1.0e2;
	841	if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
	842	if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
	843	if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
	844	if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
	845
	846	delta = min( deltamax, delta );
	847	}
	848	else
	849	{
[1126]	850	for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
[3]	851
	852	if( ts > 0.8 ) delta *= 1.e-6;
	853	if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
	854	if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
	855	if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
	856	if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
	857	if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
	858	}
	859	ts = 0.0e0;
	860	/*
	861	out = fopen( outlog, "a" );
	862	fprintf(out, " alt:%e; delta:%e; time:%e; fin:%e\n",(RA[j]-R0),delta,time,(*FIN));
	863	fclose( out );
	864	*/
	865	}
	866
	867	/* +++++++++++++++++++ */
	868	/* end of time loop. */
	869	/* +++++++++++++++++++ */
	870
	871	for( i = 0; i <= ST-1; i++ )
	872	if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
	873	fluxCH4 = ( MASS[i]/(6.022e23time) );
	874
[1126]	875	} /* boucle j */
[3]	876
	877
[1126]	878	/*
	879	==========================================================================
[3]	880
[1126]	881	FINALISATION:
	882	===============
	883	*/
	884	for( i = 0; i <= ST-1; i++ )
	885	if( strcmp(corps[i],"CH4") == 0 )
	886	fluxCH4 = ( MASS[i]/(6.022e23time) );
	887
[3]	888	/* Niveau de N2 */
	889	/* ------------ */
	890
	891	for( j = 0; j <= NLEV-1; j++ )
	892	{
	893	conv = 0.0e0;
[1126]	894	for( i = 0; i <= ST-1; i++ )
	895	if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
	896	for( i = 0; i <= ST-1; i++ )
	897	if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
[3]	898	}
	899
	900	if( (*aerprod) == 1 )
	901	{
	902	fdm2d( k_dep, 1, 5, 1 );
	903	fdm2d( faer, 1, 5, 1 );
	904	fdm1d( productaer, 0 );
	905	fdm1d( mmolaer, 0 );
	906	fdm1d( csurn, 0 );
	907	fdm1d( csurh, 0 );
	908	}
	909
	910	fdm1d( fl, 0 );
	911	fdm1d( fp, 0 );
	912	fdm1d( mu, 0 );
[1126]	913	fdm2d( ym1, 0, NC-1, 0 );
	914	fdm2d( f, 0, NC-1, 0 );
[3]	915	fdm2d( jac, 0, NC-1, 0 );
[1126]	916	fdm2d( c, 0, NLEV-1, 0 );
	917	fdm3d( a, 0, NLEV-1, 0, NC-1, 0 );
	918	fdm3d( b, 0, NLEV-1, 0, NC-1, 1 );
[3]	919	}

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

Download in other formats: