Context Navigation

gptitan.c @ 3094

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

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

Context Navigation

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

Download in other formats: