Context Navigation

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

Last change on this file since 3557 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

Line
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
13	void gptitan_(
14	double RA, double TEMP, double *NB,
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],
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];
25	int i,j,k,l;
26	int ireac,ncom1,ncom2;
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;
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 */
55	printf("CHIMIE: lat=%d\n",(*LAT)+1);
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" );
66	out = fopen( outlog, "w" );
67	fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
68	fclose( out );
69
70	deltamax = 1.e5;
71	test = 1.0e-15;
72
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
79	/* DEBUG
80	out = fopen( outlog, "a" );
81	fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
82	fclose( out );
83	*/
84	fl = dm1d( 0, NC-1 );
85	fp = dm1d( 0, NC-1 );
86	mu = dm1d( 0, NLEV-1 );
87	ym1 = dm2d( 0, NC-1, 0, NLEV-1 );
88	f = dm2d( 0, NC-1, 0, NLEV-1 );
89	jac = dm2d( 0, NC-1, 0, NC-1 );
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 );
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
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
134	/* ****************** */
135	/* Time loop: */
136	/* ****************** */
137
138	time = ts = 0.0e0;
139	delta = 1.e-3;
140
141	while( time < (*FIN) )
142	{
143
144
145	/* DEBUG
146	for( j = NLEV-1; j >= NLD; j-- )
147	{
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 );
157	}
158	exit(0);
159	*/
160
161
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]);
257	*/
258
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	{
589	time = ts = 0.0e0;
590	delta = 1.e-3;
591
592	/* ++++++++++++ */
593	/* time loop. */
594	/* ++++++++++++ */
595
596	while( time < (*FIN) )
597	{
598
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	}
606
607	/* Chimie */
608	/* ------ */
609
610	/* productions et pertes chimiques */
611	for( i = 0; i <= ST-1; i++ )
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
651
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	}
708
709
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
725	fp[utilaer[0]] += ( dyh * NB[j] );
726	/* pourquoi pas 2 ?? cf gptit dans 2da... /
727
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] );
731	/* pourquoi pas 2 ?? cf gptit dans 2da... /
732	}
733
734
735	/* Backup jacobian level j. */
736	/* ------------------------ */
737	for( i = 0; i <= ST-1; i++ )
738	{
739	for( k = 0; k <= ST-1; k++ )
740	a[j][i][k] = jac[i][k];
741	f[i][j] = fp[i] - fl[i] * Y[i][j];
742	}
743
744
745	/* finition pour inversion */
746	/* ----------------------- */
747
748	for( i = 0; i <= ST-1; i++ )
749	{
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;
756	}
757
758
759	/* Inversion of matrix cf method LU */
760	/* -------------------------------- */
761
762	/* inversion by blocs: */
763	/* Hydrocarbons */
764
765	solve_b( a, f, j, 0, NHC-1 );
766	for( i = 0; i <= NHC-1; i++ )
767	{
768	Y[i][j] += f[i][j];
769	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
770	}
771
772	/* Nitriles */
773
774	solve_b( a, f, j, NHC, ST-1 );
775	for( i = NHC+1; i <= ST-1; i++ )
776	{
777	Y[i][j] += f[i][j];
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	/* ---------- */
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	}
791
792	/* ------------------ */
793	/* Tests et evolution */
794	/* ------------------ */
795
796	/* Calcul deviation */
797	/* ---------------- */
798
799	for( i = 0; i <= ST-1; i++ )
800	{
801	test = 1.0e-15;
802	if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
803	{
804	conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
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
822	/* test deviation */
823	/* -------------- */
824
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",
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] );
834	fclose( out );
835	// exit(0);
836	Y[i][j] = 1.e-20;
837	}
838	for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
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	{
850	for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
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
875	} /* boucle j */
876
877
878	/*
879	==========================================================================
880
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
888	/* Niveau de N2 */
889	/* ------------ */
890
891	for( j = 0; j <= NLEV-1; j++ )
892	{
893	conv = 0.0e0;
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;
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 );
913	fdm2d( ym1, 0, NC-1, 0 );
914	fdm2d( f, 0, NC-1, 0 );
915	fdm2d( jac, 0, NC-1, 0 );
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 );
919	}

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

Download in other formats: