Context Navigation

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

Last change on this file since 1957 was 1956, checked in by jvatant, 7 years ago
+ Add surface methane CH4 flux pseudo-evap diagnostic and include it in bottom layer zdqcond as >0. + minor cosmetic changes in gptitan.c --JVO
File size: 28.3 KB

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