source: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ytop @ 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.9 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#include "titan.h"
10
11void 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
237printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
238if(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]);
241if(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]);
244if(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]);
247if(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]);
250if(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]);
253if(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
286printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
287if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
288printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
289if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
290printf("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.022e23*time) );
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.