source: trunk/libf/chimtitan/gptitan.c @ 16

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