source: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ftopvar @ 1058

Last change on this file since 1058 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: 15.6 KB
Line 
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     float *RA, float *TEMP, float *NB,
15     char CORPS[][10], float Y[][NLEV], float *FTOP,
16     float *DECLIN, float *FIN, int *LAT, float *MASS,
17     float *botCH4,
18     float KRPD[][NLEV][RDISS+1][15], float 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     float MAER[][NLEV], float PRODAER[][NLEV],
22     float CSN[][NLEV], float CSH[][NLEV],
23     int *htoh2, float *surfhaze)
24{
25   char   outlog[100],corps[100][10];
26   int    dec,declinint,i,j,k,l;
27   int    ireac,ncom1,ncom2;
28   float  *fl,*fp,*mu,**jac,*ym1;
29   float  *fluxtop,fluxCH4;
30   float  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
31   float  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   float  dyh,dyh2;
38
39/* va avec aer */
40   float  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
41   float  **k_dep,**faer;
42   float  *productaer,*csurn,*csurh,*mmolaer;
43
44   if( (*aerprod) == 1 )
45   {
46    k_dep = rm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
47    faer  = rm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
48    productaer  = rm1d( 0, 3 );     /* local production rate by pathways */
49    mmolaer     = rm1d( 0, 3 );     /* local molar mass by pathways */
50    csurn = rm1d( 0, 3 );           /* local C/N by pathways */
51    csurh = rm1d( 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     = rm1d( 0,   NC-1 );
74   fl      = rm1d( 0,   NC-1 );
75   fp      = rm1d( 0,   NC-1 );
76   fd      = dm1d( 0,   NC-1 );
77   mu      = rm1d( 0, NLEV-1 );
78   fluxtop = rm1d( 0,   NC-1 );
79   jac     = rm2d( 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", "w" );
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
164      time     = ts = 0.0e0;
165/*      delta    = (*FIN);  */
166      delta    = 1.e-3;
167
168      for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
169
170/* ++++++++++++ */
171/*  time loop.  */
172/* ++++++++++++ */
173
174      while( time < (*FIN) )     
175      {
176
177/* Calcul de f et de la matrice jacobienne */
178/* --------------------------------------- */
179
180         for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
181         {
182            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
183
184            fp[i] = fl[i] = 0.0e0;                    /* init for next step */
185            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
186
187            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
188            {
189               ireac = prod[i][l];                  /* Number of the reaction involves. */
190               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
191               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
192               {
193                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
194                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
195               }
196               else                                 /* General case. */
197               {
198                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
199                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
200                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
201                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
202               }
203            }
204           
205            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
206            {
207               ireac = perte[i][l][0];              /* Reaction number. */
208               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
209               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
210               {
211                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
212                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
213               }
214               else                                 /* General case. */
215               {
216                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
217                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
218                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
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
239printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
240if(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]);
243if(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]);
246if(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]);
249if(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]);
252if(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]);
255if(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/* H -> H2 on haze particles */
281/* ------------------------- */
282         if( (*htoh2) == 1 )
283         {
284              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
285/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
286
287/* DEBUG
288printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
289if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
290printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
291if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
292printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
293*/
294
295              fp[utilaer[0]] += ( dyh  * NB[j] );
296              fp[utilaer[1]] += ( dyh2 * NB[j] );
297              if( Y[utilaer[0]][j] != 0.0 )
298       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
299         }
300
301         for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
302         {
303            for( k = 0; k <= ST-1; k++ )
304            {
305               jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
306               if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
307               jacd[i][k] = (double)jac[i][k];
308            }
309           
310            fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
311         }
312
313/*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
314         {
315            Y[i][j] = ??? ;
316         }
317*/
318
319/* Inversion of matrix cf method LU */
320/* -------------------------------- */
321
322/* inversion by blocs: */
323/* Hydrocarbons */
324
325         solve( jacd, fd, 0, NHC-1 );             
326
327         for( i = 0; i <= NHC-1; i++ )
328         {
329            Y[i][j] += (float)fd[i];
330            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
331         }
332
333/* Nitriles */
334
335         solve( jacd, fd, NHC, ST-1 );             
336
337         for( i = NHC+1; i <= ST-1; i++ )
338         {
339            Y[i][j] += (float)fd[i];
340            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
341         }
342
343/* end inversion by blocs: */
344
345         for( i = 0; i <= ST-1; i++ )
346         {
347
348/* CH4 au sol */
349/* ---------- */
350
351            if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
352            {
353                fluxCH4 += (*botCH4 - Y[i][0]);
354                Y[i][0] = *botCH4;
355            }
356
357         }
358         
359/* test evolution delta */
360/* -------------------- */
361
362         for( i = 0; i <= ST-1; i++ )
363         {
364            test = 1.0e-15;
365            if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
366            {
367               conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
368
369               if( conv > ts )
370               {
371/*
372                  if( conv >= 0.1 )
373                  {
374                     out = fopen( outlog, "a" );
375                     fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
376                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
377                     fclose( out );
378                  }
379*/
380                  ts = conv;
381               }
382            }
383         }
384
385         if( ts < 0.1e0 )
386         {
387            for( i = 0; i <= ST-1; i++ )
388                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
389                 {
390                  out = fopen( outlog, "a" );
391                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
392                           corps[i], ym1[i], Y[i][j], j );
393                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
394                  fclose( out );
395                /*  exit(0); */
396                  Y[i][j] = 1.e-20;
397                 }
398            for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
399            time += delta;
400            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
401            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
402            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
403            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
404            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
405         
406            delta = min( deltamax, delta );
407         }
408         else
409         {
410            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
411
412            if(   ts > 0.8 )                    delta *= 1.e-6;
413            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
414            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
415            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
416            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
417            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
418         }
419         ts = 0.0e0;
420      }               
421
422/* +++++++++++++++++++ */
423/*  end of time loop.  */
424/* +++++++++++++++++++ */
425
426      for( i = 0; i <= ST-1; i++ )
427         if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
428            fluxCH4 *= ( MASS[i]/(6.022e23*time) );
429
430   }
431   
432/* **************** */       
433/* end of main loop */
434/* **************** */       
435
436/* Plafond: !! OU !!  flux vertical     */
437/* ------------------------------------ */
438
439   for( i = 0; i <= ST-1; i++ )
440      if( FTOP[i] != 0.0e0 )
441      {
442             fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
443             Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
444             Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
445// on ajuste aussi le niveau dans la derniere couche...
446// pour eviter les effets vers le haut
447             Y[i][NLEV-1]  = Y[i][NLEV-2];
448      }
449
450/* Niveau de N2 */
451/* ------------ */
452     
453   for( j = 0; j <= NLEV-1; j++ )
454   {
455      conv = 0.0e0;
456      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
457      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
458   }
459
460   if( (*aerprod) == 1 )
461   {
462    frm2d( k_dep, 1, 5, 1 );
463    frm2d(  faer, 1, 5, 1 );
464    frm1d( productaer,  0 );
465    frm1d( mmolaer,  0 );
466    frm1d( csurn, 0 );
467    frm1d( csurh, 0 );
468   }
469
470   frm1d(     ym1, 0 );
471   frm1d(      fl, 0 );
472   frm1d(      fp, 0 );
473   fdm1d(      fd, 0 );
474   frm1d(      mu, 0 );
475   frm1d( fluxtop, 0 );
476   frm2d(     jac, 0,   NC-1, 0 );
477   fdm2d(    jacd, 0,   NC-1, 0 );
478}
Note: See TracBrowser for help on using the repository browser.