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

Rev	Line
[3]	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
[1057]	10	void gptitan_(
[1956]	11	double RA, double TEMP, double *NB,
	12	char CORPS[][10], double Y[][NLEV],
	13	double FIN, double LAT, double *MASS, double MD[][NLEV],
[1958]	14	double *KEDD, double KRATE[][NLEV],
[1956]	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)
[3]	20	{
[2099]	21	char outlog[100];
[1956]	22	int i,j,k,l;
	23	int ireac,ncom1,ncom2;
[1958]	24	int in2, ih, ih2, in4s;
[1956]	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;
[3]	32
[1956]	33	/* va avec htoh2 */
	34	double dyh,dyh2;
[3]	35
[1956]	36	/* va avec aer */
	37	double dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
	38	double k_dep,faer;
	39	double productaer,csurn,csurh,mmolaer;
[3]	40
[1956]	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	}
[3]	50
[1956]	51	/* DEBUG */
	52	printf("CHIMIE: lat=%g\n",(*LAT));
	53	/**/
[3]	54
[1965]	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	/**/
[3]	60
[1956]	61	strcpy( outlog, "chimietitan" );
	62	strcat( outlog, ".log" );
	63	out = fopen( outlog, "w" );
	64	fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
	65	fclose( out );
[1126]	66
[1956]	67	deltamax = 1.e5;
	68	test = 1.0e-15;
[3]	69
[1956]	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;
[1126]	75
[1956]	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 );
[3]	90
[1956]	91	/* DEBUG */
	92	/*
	93	out = fopen( "err.log", "a" );
	94	fprintf( out,"%s\n", );
	95	fclose( out );
	96	*/
[3]	97
[1956]	98	/* init indexes of special species */
	99	for( i = 0; i <= ST-1; i++ )
	100	{
[2099]	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;
[1956]	105	}
	106
[1958]	107	/* initialisation mu */
[3]	108
[1956]	109	for( j = 0; j <= NLEV-1; j++ )
	110	{
[3]	111	mu[j] = 0.0e0;
	112	for( i = 0; i <= ST-1; i++ )
[1956]	113	{
	114	mu[j] += ( MASS[i] * Y[i][j] );
	115	}
	116	}
[3]	117
[1956]	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);
[1126]	121
[1956]	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	==========================================================================
[1126]	128
[1956]	129	PREMIERE ETAPE:
	130	===============
	131	INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
	132	===============
	133	*/
[1126]	134
[1956]	135	/* ****************** */
	136	/* Time loop: */
	137	/* ****************** */
[3]	138
[1956]	139	time = ts = 0.0e0;
	140	delta = 1.e-3;
[3]	141
[1956]	142	while( time < (*FIN) )
	143	{
[1126]	144
	145
[1956]	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 );
[3]	152
[1956]	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]);
[2099]	156	for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",CORPS[i],Y[i][j]);
[1956]	157	fclose( out );
	158	}
	159	exit(0);
	160	*/
[3]	161
	162
[1956]	163	/* ------------------------------ */
	164	/* Calculs variations et jacobien */
	165	/* ------------------------------ */
[1126]	166
[1956]	167	for( j = NLEV-1; j >= NLD; j-- )
	168	{
[1126]	169
[1956]	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	}
[1126]	177
[1956]	178	/* Chimie */
	179	/* ------ */
[1126]	180
[1956]	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 */
[1126]	185
[1956]	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	}
[1126]	203
[1956]	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	}
[1126]	221
	222
[1956]	223	/* Aerosols */
	224	/* -------- */
	225	if( (*aerprod) == 1 )
	226	{
[2099]	227	aer(CORPS,TEMP,NB,Y,&j,k_dep,faer,
[1956]	228	&dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
	229	mmolaer,productaer,csurn,csurh);
[1126]	230
[1956]	231	for( i = 0; i <= 3; i++ )
	232	{
	233	PRODAER[i][j] = productaer[i];
[1126]	234	MAER[i][j] = mmolaer[i];
[1956]	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.))
[2099]	241	printf("fp(%s) =%e; dyc2h2 =%e\n",CORPS[utilaer[2]],
[1956]	242	fp[utilaer[2]],dyc2h2*NB[j]);
	243	if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
[2099]	244	printf("fp(%s) =%e; dyhcn =%e\n",CORPS[utilaer[5]],
[1956]	245	fp[utilaer[5]],dyhcn*NB[j]);
	246	if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
[2099]	247	printf("fp(%s) =%e; dyhc3n =%e\n",CORPS[utilaer[6]],
[1956]	248	fp[utilaer[6]],dyhc3n*NB[j]);
	249	if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
[2099]	250	printf("fp(%s) =%e; dynccn =%e\n",CORPS[utilaer[13]],
[1956]	251	fp[utilaer[13]],dynccn*NB[j]);
	252	if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
[2099]	253	printf("fp(%s) =%e; dych3cn=%e\n",CORPS[utilaer[14]],
[1956]	254	fp[utilaer[14]],dych3cn*NB[j]);
	255	if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
[2099]	256	printf("fp(%s) =%e; dyc2h3cn=%e\n",CORPS[utilaer[15]],
[1956]	257	fp[utilaer[15]],dyc2h3cn*NB[j]);
	258	*/
[3]	259
[1956]	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	}
[1126]	279
	280
[1956]	281	/* H -> H2 on haze particles */
	282	/* ------------------------- */
	283	if( (*htoh2) == 1 )
	284	{
[2099]	285	heterohtoh2(CORPS,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
[1956]	286	/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
[1126]	287
[1956]	288	/* DEBUG
	289	printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
	290	if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
[2099]	291	printf("fp(%s) = %e; dyh = %e\n",CORPS[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
[1956]	292	if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
[2099]	293	printf("fp(%s) = %e; dyh2 = %e\n",CORPS[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
[1956]	294	*/
[1126]	295
	296	fp[utilaer[0]] += ( dyh * NB[j] );
[1956]	297	/* pourquoi pas 2 ?? cf gptit dans 2da... /
[1126]	298
	299	fp[utilaer[1]] += ( dyh2 * NB[j] );
	300	if( Y[utilaer[0]][j] != 0.0 )
[1956]	301	jac[utilaer[0]][utilaer[0]] += ( dyh * NB[j] / Y[utilaer[0]][j] );
	302	/* pourquoi pas 2 ?? cf gptit dans 2da... /
	303	}
[1126]	304
	305
[1956]	306	/* Backup jacobian level j. */
	307	/* ------------------------ */
	308	for( i = 0; i <= ST-1; i++ )
[1126]	309	for( k = 0; k <= ST-1; k++ )
[1956]	310	a[j][i][k] = jac[i][k];
[1126]	311
	312
[1956]	313	/* Diffusion verticale et flux exterieurs */
	314	/* -------------------------------------- */
[1126]	315
[1956]	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	{
[1126]	322
[1956]	323	/* First level. */
	324	if( j == NLD )
	325	{
	326	v = dv = 0.0e0;
	327	dra = RA[j+1]-RA[j];
[1126]	328
[1956]	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];
[1126]	355
[1956]	356	/* Jeans escape */
	357	if( i == ih )
	358	{
	359	dv = top_H * NB[NLEV-1]
[1126]	360	* (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
[1956]	361	v = dv * Y[i][NLEV-1];
	362	}
	363	else if( i == ih2 )
	364	{
	365	dv = top_H2 * NB[NLEV-1]
[1126]	366	* (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
[1956]	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.));
[1126]	373
[1956]	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;
[1126]	404
[1956]	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	}
[1126]	444
	445
	446
[1956]	447	/* finition pour inversion */
	448	/* ----------------------- */
[1126]	449
[1956]	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	}
[1126]	459
[1956]	460	}
[1126]	461
	462
[1956]	463	/* -------------------------------- */
	464	/* Inversion of matrix cf method LU */
	465	/* -------------------------------- */
[1126]	466
[1956]	467	for( j = NLD+1; j <= NLEV-1; j++ )
	468	{
[1965]	469	#ifdef LAPACK
	470	solve_lapack( a, j-1, 0, ST-1 );
	471	#else
[1956]	472	solve( a, j-1, 0, ST-1 );
[1965]	473	#endif
[1956]	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	}
[1965]	485	#ifdef LAPACK
	486	solve_lapack( a, NLEV-1, 0, ST-1 );
	487	#else
[1956]	488	solve( a, NLEV-1, 0, ST-1 );
[1965]	489	#endif
[1956]	490	for( j = NLEV-1; j >= NLD; j-- )
	491	{
	492	if( j != NLEV-1 )
[1126]	493	for( i = 0; i <= ST-1; i++ ) f[i][j] -= ( c[j][i] * b[j+1][i][1] );
[1956]	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	}
[1126]	503
[1956]	504	/* ------------------ */
	505	/* Tests et evolution */
	506	/* ------------------ */
[1126]	507
[1956]	508	/* Calcul deviation */
	509	/* ---------------- */
[1126]	510
[1956]	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));
[2099]	523	fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),CORPS[i],ym1[i],Y[i][j],time,delta);
[1956]	524	fclose( out );
	525	}
	526	*/
[1126]	527	ts = conv;
[1956]	528	}
	529	}
[1126]	530
[1956]	531	/* test deviation */
	532	/* -------------- */
[1126]	533
[1956]	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	{
[1126]	540	out = fopen( outlog, "a" );
	541	fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
[2099]	542	CORPS[i], ym1[i], Y[i][j], j );
[1126]	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);
[1956]	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;
[1126]	558
[1956]	559	// if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
	560	// if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
[1126]	561
[1956]	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];
[1126]	568
[1956]	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;
[1126]	577
[1956]	578	out = fopen( outlog, "a" );
	579	fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
	580	fclose( out );
[1126]	581
[1956]	582	}
	583	/* **************** */
	584	/* end of time loop */
	585	/* **************** */
[1126]	586
[1956]	587	/*
	588	==========================================================================
[1126]	589
[1956]	590	SECONDE ETAPE:
	591	===============
[1126]	592	INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
[1956]	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;
[3]	600
[1956]	601	/* ++++++++++++ */
	602	/* time loop. */
	603	/* ++++++++++++ */
[3]	604
[1956]	605	while( time < (*FIN) )
	606	{
[3]	607
[1956]	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	}
[3]	615
[1956]	616	/* Chimie */
	617	/* ------ */
[1126]	618
[1956]	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 */
[3]	623
[1956]	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	}
[3]	641
[1956]	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	}
[3]	659
[1126]	660
[1956]	661	/* Aerosols */
	662	/* -------- */
	663	if( (*aerprod) == 1 )
	664	{
[2099]	665	aer(CORPS,TEMP,NB,Y,&j,k_dep,faer,
[1956]	666	&dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
	667	mmolaer,productaer,csurn,csurh);
[3]	668
[1956]	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.))
[2099]	679	printf("fp(%s) =%e; dyc2h2 =%e\n",CORPS[utilaer[2]],
[1956]	680	fp[utilaer[2]],dyc2h2*NB[j]);
	681	if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
[2099]	682	printf("fp(%s) =%e; dyhcn =%e\n",CORPS[utilaer[5]],
[1956]	683	fp[utilaer[5]],dyhcn*NB[j]);
	684	if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
[2099]	685	printf("fp(%s) =%e; dyhc3n =%e\n",CORPS[utilaer[6]],
[1956]	686	fp[utilaer[6]],dyhc3n*NB[j]);
	687	if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
[2099]	688	printf("fp(%s) =%e; dynccn =%e\n",CORPS[utilaer[13]],
[1956]	689	fp[utilaer[13]],dynccn*NB[j]);
	690	if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
[2099]	691	printf("fp(%s) =%e; dych3cn=%e\n",CORPS[utilaer[14]],
[1956]	692	fp[utilaer[14]],dych3cn*NB[j]);
	693	if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
[2099]	694	printf("fp(%s) =%e; dyc2h3cn=%e\n",CORPS[utilaer[15]],
[1956]	695	fp[utilaer[15]],dyc2h3cn*NB[j]);
	696	*/
[3]	697
[1956]	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	}
[1126]	717
	718
[1956]	719	/* H -> H2 on haze particles */
	720	/* ------------------------- */
	721	if( (*htoh2) == 1 )
	722	{
[2099]	723	heterohtoh2(CORPS,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
[1956]	724	/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
[3]	725
[1956]	726	/* DEBUG
	727	printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
	728	if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
[2099]	729	printf("fp(%s) = %e; dyh = %e\n",CORPS[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
[1956]	730	if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
[2099]	731	printf("fp(%s) = %e; dyh2 = %e\n",CORPS[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
[1956]	732	*/
[3]	733
[1956]	734	fp[utilaer[0]] += ( dyh * NB[j] );
	735	/* pourquoi pas 2 ?? cf gptit dans 2da... /
[1126]	736
[1956]	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	}
[3]	742
[1126]	743
[1956]	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	}
[3]	752
[1126]	753
[1956]	754	/* finition pour inversion */
	755	/* ----------------------- */
[1126]	756
[1956]	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	}
[3]	766
[1126]	767
[1956]	768	/* Inversion of matrix cf method LU */
	769	/* -------------------------------- */
[3]	770
[1956]	771	/* inversion by blocs: */
	772	/* Hydrocarbons */
[3]	773
[1956]	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	}
[3]	780
[1956]	781	/* Nitriles */
[3]	782
[1956]	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	}
[3]	789
[1956]	790	/* end inversion by blocs: */
[3]	791
[1956]	792	/* ------------------ */
	793	/* Tests et evolution */
	794	/* ------------------ */
[3]	795
[1956]	796	/* Calcul deviation */
	797	/* ---------------- */
[3]	798
[1956]	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];
[3]	805
[1956]	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) );
[2099]	813	fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),CORPS[i],ym1[i],Y[i][j],time,delta);
[1956]	814	fclose( out );
	815	}
	816	*/
	817	ts = conv;
	818	}
	819	}
	820	}
[3]	821
[1956]	822	/* test deviation */
	823	/* -------------- */
[1126]	824
[1956]	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",
[2099]	832	CORPS[i], ym1[i][j], Y[i][j], j );
[1956]	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;
[3]	845
[1956]	846	delta = min( deltamax, delta );
	847	}
	848	else
	849	{
	850	for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
[3]	851
[1956]	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	}
[3]	866
[1956]	867	/* +++++++++++++++++++ */
	868	/* end of time loop. */
	869	/* +++++++++++++++++++ */
[3]	870
[1956]	871	} /* boucle j */
[3]	872
	873
[1956]	874	/*
	875	==========================================================================
[3]	876
[1956]	877	FINALISATION:
	878	===============
[3]	879
[1956]	880	/* Niveau de N2 */
	881	/* ------------ */
[1126]	882
[1956]	883	for( j = 0; j <= NLEV-1; j++ )
	884	{
[3]	885	conv = 0.0e0;
[1126]	886	for( i = 0; i <= ST-1; i++ )
[1956]	887	if( i != in2 ) conv += Y[i][j];
	888	Y[in2][j] = 1. - conv;
	889	}
[3]	890
[1956]	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	}
[3]	900
[1956]	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 );
[3]	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: