Context Navigation

gptitan.ytop @ 86

Last change on this file since 86 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

chantiers : pour communiquer sur nos projets de modifs
documentation : pour stocker les docs

Ajout de:

libf/phytitan : physique de Titan
libf/chimtitan: chimie de Titan
libf/phyvenus : physique de Venus

File size: 15.9 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	#include "titan.h"
	10
	11	void gptitan_(const int *NLAT,
	12	float RA, float TEMP, float *NB,
	13	char CORPS[][10], float Y[][NLEV], float *FTOP,
	14	float YTOP, float DECLIN, float FIN, int LAT, float *MASS,
	15	float *botCH4,
	16	float KRPD[][NLEV][RDISS+1][15], float KRATE[][NLEV],
	17	int reactif[][5], int nom_prod, int nom_perte,
	18	int prod[][200], int perte[][200][2], int aerprod, int utilaer,
	19	float MAER[][NLEV], float PRODAER[][NLEV],
	20	float CSN[][NLEV], float CSH[][NLEV],
	21	int htoh2, float surfhaze)
	22	{
	23	char outlog[100],corps[100][10];
	24	int dec,declinint,i,j,k,l;
	25	int ireac,ncom1,ncom2;
	26	float fl,fp,mu,jac,ym1;
	27	float *fluxtop,fluxCH4;
	28	float cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
	29	float rr,np,nm,factdec,s,test,time,ts,v;
	30	double fd,*jacd;
	31	char str2[15];
	32	FILE *out;
	33
	34	/* va avec htoh2 */
	35	float dyh,dyh2;
	36
	37	/* va avec aer */
	38	float dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
	39	float k_dep,faer;
	40	float productaer,csurn,csurh,mmolaer;
	41
	42	if( (*aerprod) == 1 )
	43	{
	44	k_dep = rm2d( 1, 5, 1, 3 ); /* k en s-1, reactions d'initiation */
	45	faer = rm2d( 1, 5, 1, 3 ); /* fraction de chaque compose */
	46	productaer = rm1d( 0, 3 ); /* local production rate by pathways */
	47	mmolaer = rm1d( 0, 3 ); /* local molar mass by pathways */
	48	csurn = rm1d( 0, 3 ); /* local C/N by pathways */
	49	csurh = rm1d( 0, 3 ); /* local C/H by pathways */
	50	}
	51
	52	/* DEBUG
	53	printf("CHIMIE: lat=%d declin=%e\n",(LAT)+1,(DECLIN));
	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	deltamax = 1.e5;
	65
	66	/* DEBUG
	67	out = fopen( outlog, "a" );
	68	fprintf(out,"CHIMIE: lat=%d declin=%e\n",(LAT)+1,(DECLIN));
	69	fclose( out );
	70	*/
	71	ym1 = rm1d( 0, NC-1 );
	72	fl = rm1d( 0, NC-1 );
	73	fp = rm1d( 0, NC-1 );
	74	fd = dm1d( 0, NC-1 );
	75	mu = rm1d( 0, NLEV-1 );
	76	fluxtop = rm1d( 0, NC-1 );
	77	jac = rm2d( 0, NC-1, 0, NC-1 );
	78	jacd = dm2d( 0, NC-1, 0, NC-1 );
	79
	80	/* DEBUG */
	81	/*
	82	out = fopen( "err.log", "a" );
	83	fprintf( out,"%s\n", );
	84	fclose( out );
	85	*/
	86
	87	/* initialisation krate pour dissociations */
	88
	89	if( ( (DECLIN) 10 + 267 ) < 14. )
	90	{
	91	declinint = 0;
	92	dec = 0;
	93	}
	94	else
	95	{
	96	if( ( (DECLIN) 10 + 267 ) > 520. )
	97	{
	98	declinint = 14;
	99	dec = 534;
	100	}
	101	else
	102	{
	103	declinint = 1;
	104	dec = 27;
	105	while( ( (DECLIN) 10 + 267 ) >= (float)(dec+20) )
	106	{
	107	dec += 40;
	108	declinint++;
	109	}
	110	}
	111	}
	112	if( ( (DECLIN) >= -24. ) && ( (DECLIN) <= 24. ) )
	113	factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
	114	else
	115	factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;
	116
	117	for( i = 0; i <= RDISS; i++ ) /* RDISS correspond a la dissociation de N2 par les GCR */
	118	for( j = 0; j <= NLEV-1; j++ )
	119	{
	120	if( factdec < 0. ) KRATE[i][j] = KRPD[LAT][j][i][declinint-1] fabs(factdec)
	121	+ KRPD[LAT][j][i][declinint] ( 1 + factdec);
	122	if( factdec > 0. ) KRATE[i][j] = KRPD[LAT][j][i][declinint+1] factdec
	123	+ KRPD[LAT][j][i][declinint] ( 1 - factdec);
	124	if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
	125	}
	126
	127	/* initialisation mu, CH4 au sol */
	128
	129	for( j = 0; j <= NLEV-1; j++ )
	130	{
	131	mu[j] = 0.0e0;
	132	for( i = 0; i <= ST-1; i++ )
	133	{
	134	if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
	135	{
	136	fluxCH4 = (*botCH4 - Y[i][j]);
	137	Y[i][j] = *botCH4;
	138	}
	139	mu[j] += ( MASS[i] * Y[i][j] );
	140	}
	141	}
	142
	143	/* ****************** */
	144	/* Main loop: level */
	145	/* ****************** */
	146
	147	for( j = NLEV-1; j >= 0; j-- )
	148	{
	149
	150	/* DEBUG
	151	out = fopen( outlog, "a" );
	152	fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
	153	fclose( out );
	154
	155	out = fopen( "profils.log", "w" );
	156	fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
	157	for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
	158	for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
	159	fclose( out );
	160	*/
	161
	162	time = ts = 0.0e0;
	163	/* delta = (FIN); /
	164	delta = 1.e-3;
	165
	166	for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
	167
	168	/* ++++++++++++ */
	169	/* time loop. */
	170	/* ++++++++++++ */
	171
	172	while( time < (*FIN) )
	173	{
	174
	175	/* Calcul de f et de la matrice jacobienne */
	176	/* --------------------------------------- */
	177
	178	for( i = 0; i <= ST-1; i++ ) /* productions et pertes chimiques */
	179	{
	180	Y[i][j] = max(Y[i][j],1.e-30); /* minimum */
	181
	182	fp[i] = fl[i] = 0.0e0; /* init for next step */
	183	for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
	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	/* Aerosols */
	222	/* -------- */
	223	if( (*aerprod) == 1 )
	224	{
	225	aer(corps,TEMP,NB,Y,&j,k_dep,faer,
	226	&dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
	227	mmolaer,productaer,csurn,csurh);
	228
	229	for( i = 0; i <= 3; i++ )
	230	{
	231	PRODAER[i][j] = productaer[i];
	232	MAER[i][j] = mmolaer[i];
	233	CSN[i][j] = csurn[i];
	234	CSH[i][j] = csurh[i];
	235	}
	236	/* DEBUG
	237	printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
	238	if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
	239	printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
	240	fp[utilaer[2]],dyc2h2*NB[j]);
	241	if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
	242	printf("fp(%s) =%e; dyhcn =%e\n",corps[utilaer[5]],
	243	fp[utilaer[5]],dyhcn*NB[j]);
	244	if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
	245	printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
	246	fp[utilaer[6]],dyhc3n*NB[j]);
	247	if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
	248	printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
	249	fp[utilaer[13]],dynccn*NB[j]);
	250	if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
	251	printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
	252	fp[utilaer[14]],dych3cn*NB[j]);
	253	if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
	254	printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
	255	fp[utilaer[15]],dyc2h3cn*NB[j]);
	256	*/
	257
	258	fp[utilaer[2]] -= ( dyc2h2 * NB[j] );
	259	fp[utilaer[5]] -= ( dyhcn * NB[j] );
	260	fp[utilaer[6]] -= ( dyhc3n * NB[j] );
	261	fp[utilaer[13]]-= ( dynccn * NB[j] );
	262	fp[utilaer[14]]-= ( dych3cn * NB[j] );
	263	fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
	264	if( Y[utilaer[2]][j] != 0.0 )
	265	jac[utilaer[2]][utilaer[2]] -= ( dyc2h2 * NB[j] / Y[utilaer[2]][j] );
	266	if( Y[utilaer[5]][j] != 0.0 )
	267	jac[utilaer[5]][utilaer[5]] -= ( dyhcn * NB[j] / Y[utilaer[5]][j] );
	268	if( Y[utilaer[6]][j] != 0.0 )
	269	jac[utilaer[6]][utilaer[6]] -= ( dyhc3n * NB[j] / Y[utilaer[6]][j] );
	270	if( Y[utilaer[13]][j] != 0.0 )
	271	jac[utilaer[13]][utilaer[13]] -= ( dynccn * NB[j] / Y[utilaer[13]][j] );
	272	if( Y[utilaer[14]][j] != 0.0 )
	273	jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
	274	if( Y[utilaer[15]][j] != 0.0 )
	275	jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
	276	}
	277
	278	/* H -> H2 on haze particles */
	279	/* ------------------------- */
	280	if( (*htoh2) == 1 )
	281	{
	282	heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
	283	/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
	284
	285	/* DEBUG
	286	printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
	287	if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
	288	printf("fp(%s) = %e; dyh = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
	289	if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
	290	printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
	291	*/
	292
	293	fp[utilaer[0]] += ( dyh * NB[j] );
	294	fp[utilaer[1]] += ( dyh2 * NB[j] );
	295	if( Y[utilaer[0]][j] != 0.0 )
	296	jac[utilaer[0]][utilaer[0]] += ( dyh * NB[j] / Y[utilaer[0]][j] );
	297	}
	298
	299	for( i = 0; i <= ST-1; i++ ) /* finition pour inversion */
	300	{
	301	for( k = 0; k <= ST-1; k++ )
	302	{
	303	jac[i][k] = ( -THETA delta ); /* Correction time step. */
	304	if( k == i ) jac[k][k] += NB[j]; /* Correction diagonal. */
	305	jacd[i][k] = (double)jac[i][k];
	306	}
	307
	308	fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
	309	}
	310
	311	/* for( i = ST; i <= NC-1; i++ ) pas d'inversion (soot,prod): que faire?
	312	{
	313	Y[i][j] = ??? ;
	314	}
	315	*/
	316
	317	/* Inversion of matrix cf method LU */
	318	/* -------------------------------- */
	319
	320	/* inversion by blocs: */
	321	/* Hydrocarbons */
	322
	323	solve( jacd, fd, 0, NHC-1 );
	324
	325	for( i = 0; i <= NHC-1; i++ )
	326	{
	327	Y[i][j] += (float)fd[i];
	328	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
	329	}
	330
	331	/* Nitriles */
	332
	333	solve( jacd, fd, NHC, ST-1 );
	334
	335	for( i = NHC+1; i <= ST-1; i++ )
	336	{
	337	Y[i][j] += (float)fd[i];
	338	if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
	339	}
	340
	341	/* end inversion by blocs: */
	342
	343	for( i = 0; i <= ST-1; i++ )
	344	{
	345
	346	/* CH4 au sol */
	347	/* ---------- */
	348
	349	if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
	350	{
	351	fluxCH4 += (*botCH4 - Y[i][0]);
	352	Y[i][0] = *botCH4;
	353	}
	354
	355	}
	356
	357	/* test evolution delta */
	358	/* -------------------- */
	359
	360	for( i = 0; i <= ST-1; i++ )
	361	{
	362	test = 1.0e-15;
	363	if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
	364	{
	365	conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
	366
	367	if( conv > ts )
	368	{
	369	/*
	370	if( conv >= 0.1 )
	371	{
	372	out = fopen( outlog, "a" );
	373	fprintf( out, "Lat no %d; declin:%e;", (LAT)+1, (DECLIN) );
	374	fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
	375	fclose( out );
	376	}
	377	*/
	378	ts = conv;
	379	}
	380	}
	381	}
	382
	383	if( ts < 0.1e0 )
	384	{
	385	for( i = 0; i <= ST-1; i++ )
	386	if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
	387	{
	388	out = fopen( outlog, "a" );
	389	fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
	390	corps[i], ym1[i], Y[i][j], j );
	391	for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
	392	fclose( out );
	393	/* exit(0); */
	394	Y[i][j] = 1.e-20;
	395	}
	396	for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
	397	time += delta;
	398	if( ts < 1.00e-5 ) delta *= 1.0e2;
	399	if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
	400	if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
	401	if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
	402	if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
	403
	404	delta = min( deltamax, delta );
	405	}
	406	else
	407	{
	408	for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
	409
	410	if( ts > 0.8 ) delta *= 1.e-6;
	411	if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
	412	if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
	413	if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
	414	if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
	415	if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
	416	}
	417	ts = 0.0e0;
	418	}
	419
	420	/* +++++++++++++++++++ */
	421	/* end of time loop. */
	422	/* +++++++++++++++++++ */
	423
	424	for( i = 0; i <= ST-1; i++ )
	425	if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
	426	fluxCH4 = ( MASS[i]/(6.022e23time) );
	427
	428	}
	429
	430	/* **************** */
	431	/* end of main loop */
	432	/* **************** */
	433
	434	/* Plafond: maintient de la composition */
	435	/* ------------------------------------ */
	436	/*
	437	for( i = 0; i <= ST-1; i++ )
	438	if( YTOP[i] != 0.0e0 )
	439	{
	440	fluxtop[i] = (Y[i][NLEV-1] - YTOP[i]) * MASS[i]/(6.022e23(FIN));
	441	Y[i][NLEV-1]= YTOP[i];
	442	}
	443	*/
	444	/* Plafond: !! OU !! flux vertical */
	445	/* ------------------------------------ */
	446
	447	for( i = 0; i <= ST-1; i++ )
	448	if( FTOP[i] != 0.0e0 )
	449	{
	450	fluxtop[i] = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
	451	Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2](FIN);
	452	Y[i][NLEV-2] = max(Y[i][NLEV-2],0.0e0);
	453	// on ajuste aussi le niveau dans la derniere couche...
	454	// pour eviter les effets vers le haut
	455	Y[i][NLEV-1] = Y[i][NLEV-2];
	456	}
	457
	458	/* Niveau de N2 */
	459	/* ------------ */
	460
	461	for( j = 0; j <= NLEV-1; j++ )
	462	{
	463	conv = 0.0e0;
	464	for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
	465	for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
	466	}
	467
	468	if( (*aerprod) == 1 )
	469	{
	470	frm2d( k_dep, 1, 5, 1 );
	471	frm2d( faer, 1, 5, 1 );
	472	frm1d( productaer, 0 );
	473	frm1d( mmolaer, 0 );
	474	frm1d( csurn, 0 );
	475	frm1d( csurh, 0 );
	476	}
	477
	478	frm1d( ym1, 0 );
	479	frm1d( fl, 0 );
	480	frm1d( fp, 0 );
	481	fdm1d( fd, 0 );
	482	frm1d( mu, 0 );
	483	frm1d( fluxtop, 0 );
	484	frm2d( jac, 0, NC-1, 0 );
	485	fdm2d( jacd, 0, NC-1, 0 );
	486	}

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

Download in other formats:

Original Format