Index: trunk/LMDZ.TITAN/libf/chimtitan/comp.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/comp.c	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/chimtitan/comp.c	(revision 1126)
@@ -4,8 +4,19 @@
 #include "titan.h"
 
-void comp_(char CORPS[][10], double *MASS)
+void comp_(char CORPS[][10], double *CT, double *TEMP, 
+           double *MASS, double MD[][NLEV])
 {
-   int   i;
-   char  corps[100][10];
+   int   i,j;
+   char  corps[NC+1][10];
+
+   double m,ma,epsa,sig,siga,p;
+
+   p         = 2.976e07;   /*9 10^7 R / 8 pi */
+
+      /* WARNING BACKGROUND GAS IS N2 */
+
+   ma   = 28.0134e0;               /* mass of background gas in g */
+   siga = 3.798e0;         /* Lennard-Jones length of background gas 1/10 nm */
+   epsa = 71.4e0;       /* Lennard-Jones energy of background gas */
 
    for( i = 0; i <= NC; i++)
@@ -15,4 +26,9 @@
    }
 
+   for( i = 0; i <= NC-1; i++)
+   {
+      for( j = 0; j <= NLEV-1; j++ ) MD[i][j] = 0.0e0;
+   }
+
    for( i = 0; i <= NC-1; i++ )
    {
@@ -20,4 +36,9 @@
       {
          MASS[i] = 16.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.758e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,148.6e0) );
       }
       if( strcmp(corps[i], "H") == 0 )
@@ -28,128 +49,264 @@
       {
          MASS[i] = 2.0158e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 2.827e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,59.7e0) );
       }
       if( strcmp(corps[i], "CH") == 0 )
       {
          MASS[i] = 13.02e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( ( strcmp( corps[i], "CH2" ) == 0 ) || ( strcmp( corps[i], "CH2s" ) == 0 ) )
       {
          MASS[i] = 14.03e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CH3") == 0 )
       {
          MASS[i] = 15.03e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C") == 0 )
       {
          MASS[i] = 12.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 1.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2") == 0 )
       {
          MASS[i] = 24.02e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.2e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2H") == 0 )
       {
          MASS[i] = 25.03e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2H3") == 0 )
       {
          MASS[i] = 27.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.8e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2H4") == 0 )
       {
          MASS[i] = 28.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.163e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,224.7e0) );
       }
       if( strcmp(corps[i], "C2H2") == 0 )
       {
          MASS[i] = 26.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.033e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,231.8e0) );
       }
       if( strcmp(corps[i], "C2H5") == 0 )
       {
          MASS[i] = 29.06e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2H6") == 0 )
       {
          MASS[i] = 30.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.443e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,215.7e0) );
       }
       if( strcmp(corps[i], "C3H2") == 0 )
       {
          MASS[i] = 38.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.6e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C3H3") == 0 )
       {
          MASS[i] = 39.06e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.7e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( ( strcmp(corps[i], "CH2CCH2") == 0 ) || ( strcmp(corps[i], "CH3CCH") == 0 ) )
       {
          MASS[i] = 40.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.761e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,251.8e0) );
       }
       if( strcmp(corps[i], "C3H5") == 0 )
       {
          MASS[i] = 41.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.78e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C3H6") == 0 )
       {
          MASS[i] = 42.08e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.807e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,248.9e0) );
       }
       if( strcmp(corps[i], "C3H7") == 0 )
       {
          MASS[i] = 43.09e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
        }
       if( strcmp(corps[i], "C3H8") == 0 )
       {
          MASS[i] = 44.11e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.118e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,237.1e0) );
       }
       if( strcmp(corps[i], "C4H") == 0 )
       {
          MASS[i] = 49.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.2e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( ( strcmp(corps[i], "C4H2") == 0 )||( strcmp(corps[i], "C4H2s") == 0 ) )
       {
          MASS[i] = 50.06e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.3e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C4H3") == 0 )
       {
          MASS[i] = 51.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C4H4") == 0 )
       {
          MASS[i] = 52.08e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C4H5") == 0 )
       {
          MASS[i] = 53.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C4H6") == 0 )
       {
          MASS[i] = 54.09e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.6e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C4H10") == 0 )
       {
          MASS[i] = 58.13e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.687e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,531.4e0) );
       }
       if( strcmp(corps[i], "C6H") == 0 )
       {
          MASS[i] = 73.07e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.2e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C6H2") == 0 )
       {
          MASS[i] = 74.08e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C8H2") == 0 )
       {
          MASS[i] = 98.10e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 6.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "AC6H6" ) == 0 )
       {
          MASS[i] = 78.1136e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )   /* P. G. L. */
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( ( strcmp( corps[i], "C6H5" ) == 0 ) || ( strcmp( corps[i], "AC6H5" ) == 0 ) )
       {
          MASS[i] = 77.1136e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "C6H6" ) == 0 )
       {
          MASS[i] = 78.1136e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "N2") == 0 )
@@ -160,96 +317,230 @@
       {
          MASS[i] = 14.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 1.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "NH") == 0 )
       {
          MASS[i] = 15.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CN") == 0 )
       {
          MASS[i] = 26.02e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.2e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "HCN") == 0 )
       {
          MASS[i] = 27.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.63e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,569.1e0) );
       }
       if( strcmp(corps[i], "H2CN") == 0 )
       {
          MASS[i] = 28.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.8e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "C2N") == 0 )         /* C2N */
       {
          MASS[i] = 39.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "CHCN" ) == 0 )
       {
          MASS[i]   = 39.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "CH2CN" ) == 0 )
       {
          MASS[i]   = 40.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "CH3CN" ) == 0 )
       {
          MASS[i]   = 41.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "C2H3CN" ) == 0 )
       {
          MASS[i]   = 53.06e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "NCCN") == 0 )        /* NCCN */
       {
          MASS[i] = 52.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.361e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,348.6e0) );
       }
       if( strcmp(corps[i], "C3N") == 0 )         /* C3N */
       {
          MASS[i] = 50.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "HC3N") == 0 )        /* HC3N */
       {
          MASS[i] = 51.05e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp( corps[i], "C4N2" ) == 0 )
       {
          MASS[i]   = 76.1e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "H2O") == 0 )       
       {
          MASS[i] = 18.02e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 2.641e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) 
+           / ( CT[j] * sig * omega(TEMP[j],epsa,809.1e0) );
       }
       if( ( strcmp(corps[i], "O3P") == 0 ) || ( strcmp(corps[i], "O1D") == 0 ) )      
       {
          MASS[i] = 16.0e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 1.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "OH") == 0 )        
       {
          MASS[i] = 17.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
+      }
+      if( strcmp(corps[i], "HO2") == 0 )
+      {
+         MASS[i] = 33.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
+      }
+      if( strcmp(corps[i], "H2O2") == 0 )
+      {
+         MASS[i] = 33.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
+      }
+      if( strcmp(corps[i], "O2") == 0 )
+      {
+         MASS[i] = 32.0e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
+      }
+      if( strcmp(corps[i], "O3") == 0 )
+      {
+         MASS[i] = 32.0e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CO") == 0 )          
       {
          MASS[i] = 28.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.69e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,91.7e0) );
       }
       if( strcmp(corps[i], "HCO") == 0 )        
       {
          MASS[i] = 29.02e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CO2") == 0 )       
       {
          MASS[i] = 44.01e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.941e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) 
+           / ( CT[j] * sig * omega(TEMP[j],epsa,195.2e0) );
       }
       if( strcmp(corps[i], "CH2CO") == 0 )    
       {
          MASS[i] = 42.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CH2O") == 0 )    
       {
          MASS[i] = 30.03e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.75e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( ( strcmp(corps[i], "CH2OH") == 0 ) || ( strcmp(corps[i], "CH3O") == 0 ) )  
       {
          MASS[i] = 31.04e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.4e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
       }
       if( strcmp(corps[i], "CH3OH") == 0 )       
       {
          MASS[i] = 32.042e0;
+         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
+         sig     = 1.0e-16 * pow( ( siga + 3.626e0 ), 2.0e0 );
+         for( j = 0; j <= NLEV-1; j++ )
+            MD[i][j] = sqrt( p * TEMP[j] * m )
+           / ( CT[j] * sig * omega(TEMP[j],epsa,481.8e0) );
       }
    }
Index: trunk/LMDZ.TITAN/libf/chimtitan/disso.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/disso.c	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/chimtitan/disso.c	(revision 1126)
@@ -7,5 +7,5 @@
 #include "titan.h"
 
-void disso_( double KRPD[][NLEV][RDISS+1][15], int *NLAT )
+void disso_( double KRPD[][NLRT][RDISS+1][15], int *NLAT )
 {
    static double sH2[62] = {  /* incertain en dessous de 70 et en dessus de 85... */
@@ -179,11 +179,26 @@
    double f,**flux;
    double flact;
-   char   name[60];
+   char   name[60],dir[24];
    FILE   *fp,*out;
 
-   flux = dm2d(0,NLEV,0,14);
-
+   flux = dm2d(0,NLRT,0,14);
+
+/* lecture des flux actiniques:
+   - suppose que l'executable est dans $LMDGCM/RUN/xxx/
+   - et les moyennes dans $LMDGCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
+*/
+   strcpy( dir, "../../INPUT/PHOT" );
+   if( (*NLAT) < 10 )
+   {
+     strcat( dir, "0"  );
+     strcat( dir, (const char *)ecvt((float)(*NLAT),1,&x,&x) );
+   }
+   else 
+     strcat( dir, (const char *)ecvt((float)(*NLAT),2,&x,&x) );
+   strcat( dir, "x/Moy_"  );
+   printf( "Directories for actinic fluxes: %s \n", dir );
+ 
    for( lat = 0; lat <= (*NLAT)-1; lat++ )   /* Old array is set equal to 0. */
-     for( j = 0; j <= NLEV-1; j++ )
+     for( j = 0; j <= NLRT-1; j++ )
        for( i = 0; i <= RDISS; i++ )                  
          for( s = 0; s <= 14; s++ )
@@ -195,15 +210,8 @@
    for( i = 0; i <= 16; i++ ) sC2H6[i] = 0.0e0;
 
-//   for( lat = 25; lat <= 25; lat++ )     /* Main loop on latitude */
    for( lat = 0; lat <= (*NLAT)-1; lat++ ) /*     Main loop on latitude */
    for( i = 10; i <= 310; i += 5 )             /* Main loop on wavelength. */
    {
-/* lecture des flux actiniques:
-   - suppose que l'executable est dans $LMDGCM/RUN/xxx/
-   - et les moyennes dans GCCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
-   (sauf a l'idris...)
-*/
-      strcpy( name, "../../INPUT/PHOT49/Moy_" ); 
-//    strcpy( name, "PHOT49/Moy_" ); 
+      strcpy( name, dir ); 
       if( (lat+1) < 10 )
       {
@@ -213,6 +221,6 @@
       else 
        strcat( name, (const char *)ecvt((float)(lat+1),2,&x,&x) );
-      if( i < 160 ) strcat( name, "/photgcm3a" );
-      else          strcat( name, "/photgcm3a" );
+      if( i < 160 ) strcat( name, "/photmoy3a" );
+      else          strcat( name, "/photmoy3a" );
       if( i < 100 )
       {
@@ -232,5 +240,5 @@
          exit(0);
       }
-      for( j = 0; j <= NLEV-1; j++ )
+      for( j = 0; j <= NLRT-1; j++ )
       {
          fscanf( fp,"%d ",&l );
@@ -248,5 +256,5 @@
 
       for( s = 0; s <= 14; s++ )
-       for( j = 0; j <= NLEV-1; j++ )
+       for( j = 0; j <= NLRT-1; j++ )
        {
          f = flux[j][s] * sol[l] / ( 9.5e0 * 9.5e0 );   /* !! # de reac de 0 a RDISS-1 !! */
@@ -369,11 +377,11 @@
     for( s = 0; s <= 14; s++ )
     {
-     for( j = 36; j <= NLEV-1; j++ )  /* level 36 ~ 200 km */
+     for( j = 99; j <= NLRT-1; j++ )  /* level 100 = 200 km */
       KRPD[lat][j][RDISS][s] = 1.0e-16;
-     for( j = 26; j <= 35; j++ )      /* level 26 ~ 100 km */
-      KRPD[lat][j][RDISS][s] = 1.0e-17+9.e-18*(j-26);
-     for( j = 23; j <= 25; j++ )      /* level 23 ~ 70 km */
-      KRPD[lat][j][RDISS][s] = pow(10.,(-23+2*(j-23)));
-     for( j = 0; j <= 22; j++ )
+     for( j = 49; j <= 98; j++ )      /* level 50 = 100 km */
+      KRPD[lat][j][RDISS][s] = 1.0e-17+1.8e-18*(j-49);
+     for( j = 34; j <= 48; j++ )      /* level 35 = 70 km */
+      KRPD[lat][j][RDISS][s] = pow(10.,(-23+0.4*(j-34)));
+     for( j = 0; j <= 33; j++ )
       KRPD[lat][j][RDISS][s] = 0.0e0;
     }
Index: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c	(revision 1126)
@@ -13,8 +13,7 @@
 void gptitan_(
      double *RA, double *TEMP, double *NB, 
-     char CORPS[][10], double Y[][NLEV], double *FTOP,
-     double *DECLIN, double *FIN, int *LAT, double *MASS, 
-     double *botCH4, 
-     double KRPD[][NLEV][RDISS+1][15], double KRATE[][NLEV],
+     char CORPS[][10], double Y[][NLEV],
+     double *FIN, int *LAT, double *MASS, double MD[][NLEV],
+     double *KEDD, double *botCH4, double KRATE[][NLEV],
      int reactif[][5], int *nom_prod, int *nom_perte, 
      int prod[][200], int perte[][200][2], int *aerprod, int *utilaer, 
@@ -24,11 +23,12 @@
 {
    char   outlog[100],corps[100][10];
-   int    dec,declinint,i,j,k,l;
+   int    i,j,k,l;
    int    ireac,ncom1,ncom2;
-   double  *fl,*fp,*mu,**jac,*ym1;
-   double  *fluxtop,fluxCH4;
-   double  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
-   double  rr,np,nm,factdec,s,test,time,ts,v;
-   double *fd,**jacd;
+   double  ***a,***b,**c;
+   double  *fl,*fp,*mu,**jac,**ym1,**f;
+   double  fluxCH4;
+   double  conv,delta,deltamax;
+   double  cm,cp,dim,dip,dm,dp,dym,dyp,km,kp,r,dra,dram,drap;
+   double  np,nm,s,test,time,ts,v,dv;
    char   str2[15];
    FILE   *out;
@@ -53,5 +53,5 @@
 
 /* DEBUG */
-      printf("CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
+      printf("CHIMIE: lat=%d\n",(*LAT)+1);
 /**/
 
@@ -64,19 +64,31 @@
    strcpy( outlog, "chimietitan" );
    strcat( outlog, ".log" );
+   out = fopen( outlog, "w" );
+   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
+   fclose( out );
+
    deltamax = 1.e5;
+   test = 1.0e-15; 
+
+/* valeur de r:
+r = g0 R0^2 / R * 2 * 1E-3
+avec g0 en cm/s2, R0 en km, mu et mass en g
+*/
+   r        = 21.595656e0;
 
 /*  DEBUG 
             out = fopen( outlog, "a" );
-   fprintf(out,"CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
+   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
             fclose( out );
 */
-   ym1     = dm1d( 0,   NC-1 );
    fl      = dm1d( 0,   NC-1 );
    fp      = dm1d( 0,   NC-1 );
-   fd      = dm1d( 0,   NC-1 );
    mu      = dm1d( 0, NLEV-1 );
-   fluxtop = dm1d( 0,   NC-1 );
+   ym1     = dm2d( 0,   NC-1, 0, NLEV-1 );
+   f       = dm2d( 0,   NC-1, 0, NLEV-1 );
    jac     = dm2d( 0,   NC-1, 0,   NC-1 );
-   jacd    = dm2d( 0,   NC-1, 0,   NC-1 );
+   c       = dm2d( 0, NLEV-1, 0,   NC-1 );
+   a       = dm3d( 0, NLEV-1, 0,   NC-1, 0,   NC-1 );
+   b       = dm3d( 0, NLEV-1, 0,   NC-1, 1,    2 );
 
 /* DEBUG */
@@ -87,44 +99,4 @@
 */
 
-/* initialisation krate pour dissociations */
-
-   if( ( (*DECLIN) *10 + 267 ) < 14. ) 
-   {
-     declinint = 0;
-     dec = 0;
-   }
-   else 
-   {  
-       if( ( (*DECLIN) * 10 + 267 ) > 520. ) 
-       {
-          declinint = 14;
-          dec = 534;
-       }
-       else 
-       {
-          declinint = 1;
-          dec       = 27;
-          while( ( (*DECLIN) * 10 + 267 ) >= (float)(dec+20) ) 
-          {
-            dec += 40;
-            declinint++;
-          }
-       }
-   }
-   if( ( (*DECLIN) >= -24. ) && ( (*DECLIN) <= 24. ) ) 
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
-   else
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;
-
-   for( i = 0; i <= RDISS; i++ )   /* RDISS correspond a la dissociation de N2 par les GCR */
-      for( j = 0; j <= NLEV-1; j++ )
-      { 
-         if( factdec < 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint-1] * fabs(factdec) 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 + factdec);
-         if( factdec > 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint+1] * factdec 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 - factdec);
-         if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
-      }
-        
 /* initialisation mu, CH4 au sol */
      
@@ -143,12 +115,35 @@
    }
 
+/* initialisation compo avant calcul */
+   for( j = NLEV-1; j >= 0; j-- )
+      for( i = 0; i <= ST-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
+
+/* 
+==========================================================================
+   STRATEGIE:
+    INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD 
+    PUIS INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
+==========================================================================
+
+PREMIERE ETAPE: 
+===============
+    INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD 
+===============
+*/
+
 /* ****************** */
-/*  Main loop: level  */
+/*  Time loop:        */
 /* ****************** */
 
-   for( j = NLEV-1; j >= 0; j-- )
+time     = ts = 0.0e0;
+delta    = 1.e-3;
+
+while( time < (*FIN) )     
+{
+
+
+/* DEBUG 
+   for( j = NLEV-1; j >= NLD; j-- )
    {
-
-/* DEBUG 
             out = fopen( outlog, "a" );
         fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
@@ -160,34 +155,31 @@
     for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
             fclose( out );
-
-    printf("%d %e %e %e\n",declinint,(RA[j]-R0),NB[j],TEMP[j]);
-    for (i=0;i<=RDISS-1;i++) printf("%d %e\n",i,KRPD[*LAT][j][i][declinint]);
-    for (i=0;i<=ST-1;i++)    printf("%10s %e\n",corps[i],FTOP[i]);
-
+    }
    exit(0);
 */
 
-      time     = ts = 0.0e0;
-/*      delta    = (*FIN);  */
-      delta    = 1.e-3;
-
-      for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
-
-/* ++++++++++++ */
-/*  time loop.  */
-/* ++++++++++++ */
-
-      while( time < (*FIN) )     
-      {
-
-/* Calcul de f et de la matrice jacobienne */
-/* --------------------------------------- */
-
-         for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
+
+/* ------------------------------ */
+/* Calculs variations et jacobien */
+/* ------------------------------ */
+
+   for( j = NLEV-1; j >= NLD; j-- )
+   {
+
+/* init of step */
+/* ------------ */
+         for( i = 0; i <= ST-1; i++ ) 
+         {
+            fp[i] = fl[i] = 0.0e0; 
+            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
+         }
+
+/* Chimie */
+/* ------ */
+
+/* productions et pertes chimiques */
+         for( i = 0; i <= ST-1; i++ )    
          {
             Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
-
-            fp[i] = fl[i] = 0.0e0;                    /* init for next step */
-            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
 
             for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
@@ -226,4 +218,5 @@
             }
          }
+
 
 /* Aerosols */
@@ -283,5 +276,6 @@
      jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
          }
-            
+     
+       
 /* H -> H2 on haze particles */
 /* ------------------------- */
@@ -299,39 +293,478 @@
 */
 
-              fp[utilaer[0]] += ( dyh  * NB[j] );
+              fp[utilaer[0]] += ( dyh  * NB[j] ); 
+   /* pourquoi pas *2 ?? cf gptit dans 2da... */
+
               fp[utilaer[1]] += ( dyh2 * NB[j] );
               if( Y[utilaer[0]][j] != 0.0 )
        jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
-         }
-
-         for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
-         {
+   /* pourquoi pas *2 ?? cf gptit dans 2da... */
+         }
+
+
+/* Backup jacobian level j. */
+/* ------------------------ */
+         for( i = 0; i <= ST-1; i++ )
             for( k = 0; k <= ST-1; k++ )
-            {
-               jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
-               if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
-               jacd[i][k] = (double)jac[i][k];
-            }
-            
-            fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
-         }
-
-/*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
-         {
-            Y[i][j] = ??? ;
-         }
-*/
-
+               a[j][i][k] = jac[i][k];     
+
+
+/* Diffusion verticale et flux exterieurs */
+/* -------------------------------------- */
+
+/*
+pour dy/dr, dr doit etre en cm...
+pareil pour dphi/dr
+*/
+         for( i = 0; i <= ST-1; i++ )
+         {
+
+/* First level. */
+            if( j == NLD )
+            {
+               v = dv = 0.0e0;
+               dra = RA[j+1]-RA[j];
+
+               cp  = (NB[j+1]+NB[j])/2.;  /* Mean total concentration. */
+               dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
+                     pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
+               dp  = (MD[i][j]+MD[i][j+1])/2.;     /* Mean molecular diffusion. */
+               dyp = (Y[i][j+1]-Y[i][j])/(RA[j+2]-RA[j])*2.e-5; /* Delta y level +1. */
+               kp  = (KEDD[j+1]+KEDD[j])/2.;       /* Mean eddy diffusion. */
+             /* div phi. */
+               f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp ) 
+                              + kp * dyp ) 
+                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) 
+                       + fp[i] - Y[i][j]*fl[i] + v;
+             /* dphi / dy this level. */
+               a[j][i][i] += ( cp * ( dp * 0.5e0 * dip 
+                                    - 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) ) 
+                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) + dv );
+             /* dphi / dy level +1. */
+               c[j][i]     = -THETA * delta 
+                             * cp * ( dp * 0.5e0 * dip 
+                                    + 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) ) 
+                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.));
+            }
+/* Last level. */
+            else if( j == NLEV-1 )
+            {
+               v = dv = 0.0e0;
+               dra = RA[NLEV-1]-RA[NLEV-2];
+
+   /* Jeans escape */
+               if( strcmp(corps[i], "H") == 0 )  
+               {
+                 dv = top_H  * NB[NLEV-1] 
+                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.)); 
+                 v  = dv * Y[i][NLEV-1];
+               }
+               if( strcmp(corps[i], "H2") == 0 )
+               {
+                 dv = top_H2 * NB[NLEV-1]
+                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.)); 
+                 v  = dv * Y[i][NLEV-1];
+               }
+   /* Input flux for N(4S) */
+               if( strcmp(corps[i], "N4S") == 0 )
+                 v  = top_N4S
+                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.)); 
+
+               cm  = (NB[NLEV-1]+NB[NLEV-2])/2.;  /* Mean total concentration. */
+               dim = r * (MASS[i]-(mu[NLEV-1]+mu[NLEV-2])/2.)
+                       / (TEMP[NLEV-1]+TEMP[NLEV-2]) 
+                       / pow( RA[NLEV-1],   2.0e0 );  /* Delta i,j level -1. */
+               dm  = (MD[i][NLEV-1]+MD[i][NLEV-2])/2.;    /* Mean molecular diffusion. */
+               dym = (Y[i][NLEV-1]-Y[i][NLEV-2])/dra*1.e-5; /* Delta y level -1. */
+               km  = (KEDD[NLEV-1]+KEDD[NLEV-2])/2.;      /* Mean eddy diffusion. */
+             /* div phi. */
+               f[i][NLEV-1] = fp[i] - Y[i][NLEV-1]*fl[i] - v 
+                       - cm * ( dm * ( (Y[i][NLEV-1]+Y[i][NLEV-2])/2. * dim + dym ) 
+                              + km * dym ) 
+                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)); 
+             /* dphi / dy this level */
+               a[NLEV-1][i][i] -= ( cm * ( dm * 0.5e0 * dim 
+                                    + 1.e-5/dra * (dm + km ) )
+                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)) + dv );
+             /* dphi / dy level -1. */
+               b[NLEV-1][i][2]  =  THETA * delta 
+                                  * cm * ( dm * 0.5e0 * dim 
+                                    - 1.e-5/dra * (dm + km ) )
+                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
+            }
+            else
+            {
+               v = dv = 0.0e0;
+               dram=(RA[j+1]-RA[j-1])/2.;
+               if (j<NLEV-2)
+                 drap=(RA[j+1]-RA[j-1])/2.;
+               else
+                 drap=dram;
+
+               cm  = (NB[j]+NB[j-1])/2.;       /* Mean concentration level -1. */
+               cp  = (NB[j]+NB[j+1])/2.;       /* Mean concentration level +1. */
+               dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
+                     pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
+               dim = r * (MASS[i]-(mu[j]+mu[j-1])/2.) / (TEMP[j]+TEMP[j-1]) /
+                     pow( RA[j],   2.0e0 );    /* Delta i,j level -1. */
+               dm  = (MD[i][j-1]+MD[i][j])/2.;    /* Mean molecular diffusion level -1. */
+               dp  = (MD[i][j+1]+MD[i][j])/2.;    /* Mean molecular diffusion level +1. */
+               dym = (Y[i][j]-Y[i][j-1])/dram*1.e-5; /* Delta y level -1. */
+               dyp = (Y[i][j+1]-Y[i][j])/drap*1.e-5; /* Delta y level +1. */
+               km  = (KEDD[j]+KEDD[j-1])/2.;      /* Mean eddy diffusion level -1. */
+               kp  = (KEDD[j]+KEDD[j+1])/2.;      /* Mean eddy diffusion level +1. */
+             /* div phi. */
+               f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp ) 
+                              + kp * dyp ) 
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.)) 
+                       - cm * ( dm * ( (Y[i][j]+Y[i][j-1])/2. * dim + dym ) 
+                              + km * dym ) 
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) 
+                       + fp[i] - fl[i] * Y[i][j] + v;
+             /* dphi / dy this level */
+               a[j][i][i] += ( cp * ( dp * 0.5e0 * dip 
+                                    - 1.e-5/drap * (dp + kp) ) 
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
+                             - cm * ( dm * 0.5e0 * dim 
+                                    + 1.e-5/dram * (dm + km ) )
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) );
+             /* dphi / dy level -1. */
+               b[j][i][2]  =  THETA * delta 
+                             * cm * ( dm * 0.5e0 * dim
+                                    - 1.e-5/dram * (dm + km ) )
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.));
+             /* dphi / dy level +1. */
+               c[j][i]     = -THETA * delta
+                             * cp * ( dp * 0.5e0 * dip
+                                    + 1.e-5/drap * (dp + kp) ) 
+                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.));
+            }
+         }
+
+
+
+/* finition pour inversion */
+/* ----------------------- */
+
+         for( i = 0; i <= ST-1; i++ )
+         {
+            for( k = 0; k <= ST-1; k++ )
+            {
+               a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
+               if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
+            }
+            f[i][j] *= delta;
+         }
+
+   }
+
+
+/* -------------------------------- */
 /* Inversion of matrix cf method LU */
 /* -------------------------------- */
 
+   for( j = NLD+1; j <= NLEV-1; j++ )
+   {
+         solve( a, j-1, 0, ST-1 );
+         for( i = 0; i <= ST-1; i++ )
+         {
+            s = 0.0e0;
+            for( k = 0; k <= ST-1; k++ )
+            {
+               a[j][i][k] -= ( b[j][i][2] * c[j-1][k] * a[j-1][i][k] );
+               s          += ( b[j][i][2] * f[k][j-1] * a[j-1][i][k] );
+            }
+            f[i][j] -= s;
+         }
+   }
+   solve( a, NLEV-1, 0, ST-1 );
+   for( j = NLEV-1; j >= NLD; j-- )     
+   {
+         if( j != NLEV-1 )
+            for( i = 0; i <= ST-1; i++ ) f[i][j] -= ( c[j][i] * b[j+1][i][1] );
+         for( i = 0; i <= ST-1; i++ )
+         {
+            s = 0.0e0;
+            for( k = 0; k <= ST-1; k++ ) s += ( a[j][i][k] * f[k][j] );
+            b[j][i][1]  = s;
+            Y[i][j]    += s;
+            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
+         }
+   }
+
+/* ------------------ */
+/* Tests et evolution */
+/* ------------------ */
+
+/* Calcul deviation */
+/* ---------------- */
+
+   for( j = NLD; j <= NLEV-1; j++ )
+      for( i = 0; i <= ST-1; i++ )
+         if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
+         {
+               conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
+               if( conv > ts )
+               {
+/*
+                  if( conv >= 0.1 )
+                  {
+                     out = fopen( outlog, "a" );
+                     fprintf( out, "Lat no %d;", (*LAT)+1);
+                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
+                     fclose( out );
+                  }
+*/
+                  ts = conv;
+               }
+         }
+
+/* test deviation */
+/* -------------- */
+
+         if( ts < 0.1e0 )
+         {
+            for( i = 0; i <= ST-1; i++ )
+               for( j = NLD; j <= NLEV-1; j++ )
+                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
+                 {
+                  out = fopen( outlog, "a" );
+                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
+                           corps[i], ym1[i], Y[i][j], j );
+                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
+                  fclose( out );
+                  exit(0); 
+//                  Y[i][j] = 1.e-20;
+                 }
+            for( j = NLD; j <= NLEV-1; j++ )
+               for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
+            time += delta;
+            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
+            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
+            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
+            if( ( ts > 1.00e-3 ) && ( ts < 5.0e-3 ) ) delta *= 3.0e0;
+            if( ( ts > 5.00e-3 ) && ( ts < 0.01e0 ) ) delta *= 1.5e0;
+            if( ( ts > 0.010e0 ) && ( ts < 0.03e0 ) ) delta *= 1.2e0;
+            if( ( ts > 0.030e0 ) && ( ts < 0.05e0 ) ) delta *= 1.1e0;
+
+//            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
+//            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
+         
+            delta = min( deltamax, delta );
+         }
+         else
+         {
+            for( j = NLD; j <= NLEV-1; j++ )
+               for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
+
+            if(   ts > 0.8 )                    delta *= 1.e-6;
+            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
+            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
+            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
+            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
+            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
+         }
+         ts = 0.0e0;
+
+         out = fopen( outlog, "a" );
+         fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
+         fclose( out );
+
+}
+/* **************** */        
+/* end of time loop */
+/* **************** */        
+
+/*
+==========================================================================
+
+SECONDE ETAPE: 
+===============
+    INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
+===============
+*/
+   if( NLD != 0 ) 
+   for( j = NLD-1; j >= 0; j-- )
+   {
+      time     = ts = 0.0e0;
+      delta    = 1.e-3;
+
+/* ++++++++++++ */
+/*  time loop.  */
+/* ++++++++++++ */
+
+      while( time < (*FIN) )     
+      {
+
+/* init of step */
+/* ------------ */
+         for( i = 0; i <= ST-1; i++ ) 
+         {
+            fp[i] = fl[i] = 0.0e0; 
+            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
+         }
+
+/* Chimie */
+/* ------ */
+
+/* productions et pertes chimiques */
+         for( i = 0; i <= ST-1; i++ )    
+         {
+            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
+
+            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
+            {
+               ireac = prod[i][l];                  /* Number of the reaction involves. */
+               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
+               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
+               {
+                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
+                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
+               }
+               else                                 /* General case. */
+               {
+                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
+                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
+                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
+                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
+               }
+            }
+            
+            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
+            {
+               ireac = perte[i][l][0];              /* Reaction number. */
+               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
+               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
+               {
+                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
+                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
+               }
+               else                                 /* General case. */
+               {
+                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
+                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
+                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
+               }
+            }
+         }
+
+
+/* Aerosols */
+/* -------- */
+         if( (*aerprod) == 1 )
+         {
+             aer(corps,TEMP,NB,Y,&j,k_dep,faer,
+              &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
+              mmolaer,productaer,csurn,csurh);
+
+             for( i = 0; i <= 3; i++ )
+             {
+               PRODAER[i][j] = productaer[i];
+                  MAER[i][j] = mmolaer[i];
+                   CSN[i][j] = csurn[i];
+                   CSH[i][j] = csurh[i];
+             }
+/* DEBUG
+printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
+if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
+      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
+              fp[utilaer[2]],dyc2h2*NB[j]);
+if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
+      printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
+              fp[utilaer[5]],dyhcn*NB[j]);
+if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
+      printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
+              fp[utilaer[6]],dyhc3n*NB[j]);
+if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
+      printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
+              fp[utilaer[13]],dynccn*NB[j]);
+if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
+      printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
+              fp[utilaer[14]],dych3cn*NB[j]);
+if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
+      printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
+              fp[utilaer[15]],dyc2h3cn*NB[j]);
+*/
+
+             fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
+             fp[utilaer[5]] -= (    dyhcn * NB[j] );
+             fp[utilaer[6]] -= (   dyhc3n * NB[j] );
+             fp[utilaer[13]]-= (   dynccn * NB[j] );
+             fp[utilaer[14]]-= (  dych3cn * NB[j] );
+             fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
+             if( Y[utilaer[2]][j]  != 0.0 )
+       jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
+             if( Y[utilaer[5]][j]  != 0.0 )
+       jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
+             if( Y[utilaer[6]][j]  != 0.0 )
+       jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
+             if( Y[utilaer[13]][j] != 0.0 )
+     jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
+             if( Y[utilaer[14]][j] != 0.0 )
+     jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
+             if( Y[utilaer[15]][j] != 0.0 )
+     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
+         }
+     
+       
+/* H -> H2 on haze particles */
+/* ------------------------- */
+         if( (*htoh2) == 1 )
+         {
+              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
+/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
+
+/* DEBUG 
+printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
+if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
+printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
+if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
+printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
+*/
+
+              fp[utilaer[0]] += ( dyh  * NB[j] ); 
+   /* pourquoi pas *2 ?? cf gptit dans 2da... */
+
+              fp[utilaer[1]] += ( dyh2 * NB[j] );
+              if( Y[utilaer[0]][j] != 0.0 )
+       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
+   /* pourquoi pas *2 ?? cf gptit dans 2da... */
+         }
+
+
+/* Backup jacobian level j. */
+/* ------------------------ */
+         for( i = 0; i <= ST-1; i++ )
+         {
+            for( k = 0; k <= ST-1; k++ )
+               a[j][i][k] = jac[i][k];     
+            f[i][j] = fp[i] - fl[i] * Y[i][j];
+         }
+
+
+/* finition pour inversion */
+/* ----------------------- */
+
+         for( i = 0; i <= ST-1; i++ )
+         {
+            for( k = 0; k <= ST-1; k++ )
+            {
+               a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
+               if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
+            }
+            f[i][j] *= delta;
+         }
+
+
+/* Inversion of matrix cf method LU */
+/* -------------------------------- */
+
 /* inversion by blocs: */
 /* Hydrocarbons */
 
-         solve( jacd, fd, 0, NHC-1 );             
-
+         solve_b( a, f, j, 0, NHC-1 );             
          for( i = 0; i <= NHC-1; i++ )
          {
-            Y[i][j] += (float)fd[i];
+            Y[i][j] += f[i][j];
             if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
          }
@@ -339,37 +772,35 @@
 /* Nitriles */
 
-         solve( jacd, fd, NHC, ST-1 );             
-
+         solve_b( a, f, j, NHC, ST-1 );             
          for( i = NHC+1; i <= ST-1; i++ )
          {
-            Y[i][j] += (float)fd[i];
+            Y[i][j] += f[i][j];
             if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
          }
 
 /* end inversion by blocs: */
-
-         for( i = 0; i <= ST-1; i++ )
-         {
 
 /* CH4 au sol */
 /* ---------- */
-
-            if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
-            {
-                fluxCH4 += (*botCH4 - Y[i][0]);
-                Y[i][0] = *botCH4;
-            }
-
-         }
-         
-/* test evolution delta */
-/* -------------------- */
+   for( i = 0; i <= ST-1; i++ )
+     if( ( strcmp(corps[i], "CH4") == 0 ) && (j==0) && ( Y[i][0] < *botCH4 ) )
+     {
+          fluxCH4 += (*botCH4 - Y[i][0]);
+           Y[i][0] = *botCH4;
+     }
+
+/* ------------------ */
+/* Tests et evolution */
+/* ------------------ */
+
+/* Calcul deviation */
+/* ---------------- */
 
          for( i = 0; i <= ST-1; i++ )
          {
             test = 1.0e-15; 
-            if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
-            {
-               conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
+            if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
+            {
+               conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
 
                if( conv > ts )
@@ -389,4 +820,7 @@
          }
 
+/* test deviation */
+/* -------------- */
+
          if( ts < 0.1e0 )
          {
@@ -396,11 +830,11 @@
                   out = fopen( outlog, "a" );
                   fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
-                           corps[i], ym1[i], Y[i][j], j );
-                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
+                           corps[i], ym1[i][j], Y[i][j], j );
+                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i][j],Y[i][k] );
                   fclose( out );
              //     exit(0); 
                   Y[i][j] = 1.e-20;
                  }
-            for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
+            for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
             time += delta;
             if(   ts < 1.00e-5 )                      delta *= 1.0e2;
@@ -414,5 +848,5 @@
          else
          {
-            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
+            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
 
             if(   ts > 0.8 )                    delta *= 1.e-6;
@@ -439,23 +873,16 @@
             fluxCH4 *= ( MASS[i]/(6.022e23*time) ); 
 
-   }
-   
-/* **************** */        
-/* end of main loop */
-/* **************** */        
-
-/* Plafond: !! OU !!  flux vertical     */
-/* ------------------------------------ */
-
-   for( i = 0; i <= ST-1; i++ ) 
-      if( FTOP[i] != 0.0e0 )
-      {
-             fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
-             Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
-             Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
-// on ajuste aussi le niveau dans la derniere couche...
-// pour eviter les effets vers le haut
-             Y[i][NLEV-1]  = Y[i][NLEV-2];
-      }
+   }  /*  boucle j */
+
+
+/*
+==========================================================================
+
+FINALISATION: 
+===============
+*/
+      for( i = 0; i <= ST-1; i++ ) 
+         if( strcmp(corps[i],"CH4") == 0 )
+            fluxCH4 *= ( MASS[i]/(6.022e23*time) ); 
 
 /* Niveau de N2 */
@@ -465,6 +892,8 @@
    {
       conv = 0.0e0;
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
+      for( i = 0; i <= ST-1; i++ ) 
+         if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
+      for( i = 0; i <= ST-1; i++ ) 
+         if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
    }
 
@@ -479,11 +908,12 @@
    }
 
-   fdm1d(     ym1, 0 );
    fdm1d(      fl, 0 );
    fdm1d(      fp, 0 );
-   fdm1d(      fd, 0 );
    fdm1d(      mu, 0 );
-   fdm1d( fluxtop, 0 );
+   fdm2d(     ym1, 0,   NC-1, 0 );
+   fdm2d(       f, 0,   NC-1, 0 );
    fdm2d(     jac, 0,   NC-1, 0 );
-   fdm2d(    jacd, 0,   NC-1, 0 );
+   fdm2d(       c, 0, NLEV-1, 0 );
+   fdm3d(       a, 0, NLEV-1, 0,   NC-1, 0 );
+   fdm3d(       b, 0, NLEV-1, 0,   NC-1, 1 );
 }
Index: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ftopvar
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ftopvar	(revision 1121)
+++ 	(revision )
@@ -1,478 +1,0 @@
-/* gptitan: photochimie */
-/* GCCM */
-
-/* tout est passe en simple precision */
-/* sauf pour l'inversion de la matrice */
-
-/* nitriles et hydrocarbures separes pour l'inversion */
-
-/* flux variable au sommet */
-
-#include "titan.h"
-
-void gptitan_(const int *NLAT,
-     float *RA, float *TEMP, float *NB, 
-     char CORPS[][10], float Y[][NLEV], float *FTOP,
-     float *DECLIN, float *FIN, int *LAT, float *MASS, 
-     float *botCH4, 
-     float KRPD[][NLEV][RDISS+1][15], float KRATE[][NLEV],
-     int reactif[][5], int *nom_prod, int *nom_perte, 
-     int prod[][200], int perte[][200][2], int *aerprod, int *utilaer, 
-     float MAER[][NLEV], float PRODAER[][NLEV], 
-     float CSN[][NLEV], float CSH[][NLEV],
-     int *htoh2, float *surfhaze)
-{
-   char   outlog[100],corps[100][10];
-   int    dec,declinint,i,j,k,l;
-   int    ireac,ncom1,ncom2;
-   float  *fl,*fp,*mu,**jac,*ym1;
-   float  *fluxtop,fluxCH4;
-   float  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
-   float  rr,np,nm,factdec,s,test,time,ts,v;
-   double *fd,**jacd;
-   char   str2[15];
-   FILE   *out;
-
-/* va avec htoh2 */
-   float  dyh,dyh2;
-
-/* va avec aer */
-   float  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
-   float  **k_dep,**faer;
-   float  *productaer,*csurn,*csurh,*mmolaer;
-
-   if( (*aerprod) == 1 )
-   {
-    k_dep = rm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
-    faer  = rm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
-    productaer  = rm1d( 0, 3 );     /* local production rate by pathways */
-    mmolaer     = rm1d( 0, 3 );     /* local molar mass by pathways */
-    csurn = rm1d( 0, 3 );           /* local C/N by pathways */
-    csurh = rm1d( 0, 3 );           /* local C/H by pathways */
-   }
-
-/* DEBUG
-      printf("CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
-*/
-
-   for( i = 0; i <= NC; i++)
-   {
-     strcpy( corps[i], CORPS[i] );
-     corps[i][strcspn(CORPS[i], " ")] = '\0';
-   }
-   
-   strcpy( outlog, "chimietitan" );
-   strcat( outlog, ".log" );
-   deltamax = 1.e5;
-
-/*  DEBUG 
-            out = fopen( outlog, "a" );
-   fprintf(out,"CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
-            fclose( out );
-*/
-   ym1     = rm1d( 0,   NC-1 );
-   fl      = rm1d( 0,   NC-1 );
-   fp      = rm1d( 0,   NC-1 );
-   fd      = dm1d( 0,   NC-1 );
-   mu      = rm1d( 0, NLEV-1 );
-   fluxtop = rm1d( 0,   NC-1 );
-   jac     = rm2d( 0,   NC-1, 0,   NC-1 );
-   jacd    = dm2d( 0,   NC-1, 0,   NC-1 );
-
-/* DEBUG */
-/* 
-            out = fopen( "err.log", "a" );
-            fprintf( out,"%s\n", );
-            fclose( out );
-*/
-
-/* initialisation krate pour dissociations */
-
-   if( ( (*DECLIN) *10 + 267 ) < 14. ) 
-   {
-     declinint = 0;
-     dec = 0;
-   }
-   else 
-   {  
-       if( ( (*DECLIN) * 10 + 267 ) > 520. ) 
-       {
-          declinint = 14;
-          dec = 534;
-       }
-       else 
-       {
-          declinint = 1;
-          dec       = 27;
-          while( ( (*DECLIN) * 10 + 267 ) >= (float)(dec+20) ) 
-          {
-            dec += 40;
-            declinint++;
-          }
-       }
-   }
-   if( ( (*DECLIN) >= -24. ) && ( (*DECLIN) <= 24. ) ) 
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
-   else
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;
-
-   for( i = 0; i <= RDISS; i++ )   /* RDISS correspond a la dissociation de N2 par les GCR */
-      for( j = 0; j <= NLEV-1; j++ )
-      { 
-         if( factdec < 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint-1] * fabs(factdec) 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 + factdec);
-         if( factdec > 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint+1] * factdec 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 - factdec);
-         if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
-      }
-        
-/* initialisation mu, CH4 au sol */
-     
-   for( j = 0; j <= NLEV-1; j++ )
-   {
-      mu[j] = 0.0e0;
-      for( i = 0; i <= ST-1; i++ ) 
-      {
-         if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
-         {
-             fluxCH4 = (*botCH4 - Y[i][j]);
-             Y[i][j] = *botCH4;
-         }
-         mu[j] += ( MASS[i] * Y[i][j] );
-      }
-   }
-
-/* ****************** */
-/*  Main loop: level  */
-/* ****************** */
-
-   for( j = NLEV-1; j >= 0; j-- )
-   {
-
-/* DEBUG
-            out = fopen( outlog, "a" );
-        fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
-            fclose( out );
-
-            out = fopen( "profils.log", "w" );
-    fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
-    for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
-    for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
-            fclose( out );
-*/
-
-      time     = ts = 0.0e0;
-/*      delta    = (*FIN);  */
-      delta    = 1.e-3;
-
-      for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
-
-/* ++++++++++++ */
-/*  time loop.  */
-/* ++++++++++++ */
-
-      while( time < (*FIN) )     
-      {
-
-/* Calcul de f et de la matrice jacobienne */
-/* --------------------------------------- */
-
-         for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
-         {
-            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
-
-            fp[i] = fl[i] = 0.0e0;                    /* init for next step */
-            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
-
-            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
-            {
-               ireac = prod[i][l];                  /* Number of the reaction involves. */
-               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
-               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
-               {
-                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
-                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
-               }
-               else                                 /* General case. */
-               {
-                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
-                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
-                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
-                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
-               }
-            }
-            
-            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
-            {
-               ireac = perte[i][l][0];              /* Reaction number. */
-               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
-               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
-               {
-                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
-                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
-               }
-               else                                 /* General case. */
-               {
-                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
-                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
-                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
-               }
-            }
-         }
-
-/* Aerosols */
-/* -------- */
-         if( (*aerprod) == 1 )
-         {
-             aer(corps,TEMP,NB,Y,&j,k_dep,faer,
-              &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
-              mmolaer,productaer,csurn,csurh);
-
-             for( i = 0; i <= 3; i++ )
-             {
-               PRODAER[i][j] = productaer[i];
-                  MAER[i][j] = mmolaer[i];
-                   CSN[i][j] = csurn[i];
-                   CSH[i][j] = csurh[i];
-             }
-/* DEBUG
-printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
-if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
-      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
-              fp[utilaer[2]],dyc2h2*NB[j]);
-if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
-      printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
-              fp[utilaer[5]],dyhcn*NB[j]);
-if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
-      printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
-              fp[utilaer[6]],dyhc3n*NB[j]);
-if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
-      printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
-              fp[utilaer[13]],dynccn*NB[j]);
-if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
-      printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
-              fp[utilaer[14]],dych3cn*NB[j]);
-if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
-      printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
-              fp[utilaer[15]],dyc2h3cn*NB[j]);
-*/
-
-             fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
-             fp[utilaer[5]] -= (    dyhcn * NB[j] );
-             fp[utilaer[6]] -= (   dyhc3n * NB[j] );
-             fp[utilaer[13]]-= (   dynccn * NB[j] );
-             fp[utilaer[14]]-= (  dych3cn * NB[j] );
-             fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
-             if( Y[utilaer[2]][j]  != 0.0 )
-       jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
-             if( Y[utilaer[5]][j]  != 0.0 )
-       jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
-             if( Y[utilaer[6]][j]  != 0.0 )
-       jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
-             if( Y[utilaer[13]][j] != 0.0 )
-     jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
-             if( Y[utilaer[14]][j] != 0.0 )
-     jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
-             if( Y[utilaer[15]][j] != 0.0 )
-     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
-         }
-            
-/* H -> H2 on haze particles */
-/* ------------------------- */
-         if( (*htoh2) == 1 )
-         {
-              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
-/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
-
-/* DEBUG 
-printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
-if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
-printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
-if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
-printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
-*/
-
-              fp[utilaer[0]] += ( dyh  * NB[j] );
-              fp[utilaer[1]] += ( dyh2 * NB[j] );
-              if( Y[utilaer[0]][j] != 0.0 )
-       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
-         }
-
-         for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
-         {
-            for( k = 0; k <= ST-1; k++ )
-            {
-               jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
-               if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
-               jacd[i][k] = (double)jac[i][k];
-            }
-            
-            fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
-         }
-
-/*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
-         {
-            Y[i][j] = ??? ;
-         }
-*/
-
-/* Inversion of matrix cf method LU */
-/* -------------------------------- */
-
-/* inversion by blocs: */
-/* Hydrocarbons */
-
-         solve( jacd, fd, 0, NHC-1 );             
-
-         for( i = 0; i <= NHC-1; i++ )
-         {
-            Y[i][j] += (float)fd[i];
-            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
-         }
-
-/* Nitriles */
-
-         solve( jacd, fd, NHC, ST-1 );             
-
-         for( i = NHC+1; i <= ST-1; i++ )
-         {
-            Y[i][j] += (float)fd[i];
-            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
-         }
-
-/* end inversion by blocs: */
-
-         for( i = 0; i <= ST-1; i++ )
-         {
-
-/* CH4 au sol */
-/* ---------- */
-
-            if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
-            {
-                fluxCH4 += (*botCH4 - Y[i][0]);
-                Y[i][0] = *botCH4;
-            }
-
-         }
-         
-/* test evolution delta */
-/* -------------------- */
-
-         for( i = 0; i <= ST-1; i++ )
-         {
-            test = 1.0e-15; 
-            if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
-            {
-               conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
-
-               if( conv > ts )
-               {
-/*
-                  if( conv >= 0.1 )
-                  {
-                     out = fopen( outlog, "a" );
-                     fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
-                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
-                     fclose( out );
-                  }
-*/
-                  ts = conv;
-               }
-            }
-         }
-
-         if( ts < 0.1e0 )
-         {
-            for( i = 0; i <= ST-1; i++ )
-                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
-                 {
-                  out = fopen( outlog, "a" );
-                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
-                           corps[i], ym1[i], Y[i][j], j );
-                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
-                  fclose( out );
-                /*  exit(0); */
-                  Y[i][j] = 1.e-20;
-                 }
-            for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
-            time += delta;
-            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
-            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
-            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
-            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
-            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
-         
-            delta = min( deltamax, delta );
-         }
-         else
-         {
-            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
-
-            if(   ts > 0.8 )                    delta *= 1.e-6;
-            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
-            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
-            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
-            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
-            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
-         }
-         ts = 0.0e0;
-      }               
-
-/* +++++++++++++++++++ */
-/*  end of time loop.  */
-/* +++++++++++++++++++ */
-
-      for( i = 0; i <= ST-1; i++ ) 
-         if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
-            fluxCH4 *= ( MASS[i]/(6.022e23*time) ); 
-
-   }
-   
-/* **************** */        
-/* end of main loop */
-/* **************** */        
-
-/* Plafond: !! OU !!  flux vertical     */
-/* ------------------------------------ */
-
-   for( i = 0; i <= ST-1; i++ ) 
-      if( FTOP[i] != 0.0e0 )
-      {
-             fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
-             Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
-             Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
-// on ajuste aussi le niveau dans la derniere couche...
-// pour eviter les effets vers le haut
-             Y[i][NLEV-1]  = Y[i][NLEV-2];
-      }
-
-/* Niveau de N2 */
-/* ------------ */
-      
-   for( j = 0; j <= NLEV-1; j++ ) 
-   {
-      conv = 0.0e0;
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
-   }
-
-   if( (*aerprod) == 1 )
-   {
-    frm2d( k_dep, 1, 5, 1 );
-    frm2d(  faer, 1, 5, 1 );
-    frm1d( productaer,  0 );
-    frm1d( mmolaer,  0 );
-    frm1d( csurn, 0 );
-    frm1d( csurh, 0 );
-   }
-
-   frm1d(     ym1, 0 );
-   frm1d(      fl, 0 );
-   frm1d(      fp, 0 );
-   fdm1d(      fd, 0 );
-   frm1d(      mu, 0 );
-   frm1d( fluxtop, 0 );
-   frm2d(     jac, 0,   NC-1, 0 );
-   fdm2d(    jacd, 0,   NC-1, 0 );
-}
Index: trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ytop
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/gptitan.ytop	(revision 1121)
+++ 	(revision )
@@ -1,486 +1,0 @@
-/* gptitan: photochimie */
-/* GCCM */
-
-/* tout est passe en simple precision */
-/* sauf pour l'inversion de la matrice */
-
-/* nitriles et hydrocarbures separes pour l'inversion */
-
-#include "titan.h"
-
-void gptitan_(const int *NLAT,
-     float *RA, float *TEMP, float *NB, 
-     char CORPS[][10], float Y[][NLEV], float *FTOP,
-     float *YTOP, float *DECLIN, float *FIN, int *LAT, float *MASS, 
-     float *botCH4, 
-     float KRPD[][NLEV][RDISS+1][15], float KRATE[][NLEV],
-     int reactif[][5], int *nom_prod, int *nom_perte, 
-     int prod[][200], int perte[][200][2], int *aerprod, int *utilaer, 
-     float MAER[][NLEV], float PRODAER[][NLEV], 
-     float CSN[][NLEV], float CSH[][NLEV],
-     int *htoh2, float *surfhaze)
-{
-   char   outlog[100],corps[100][10];
-   int    dec,declinint,i,j,k,l;
-   int    ireac,ncom1,ncom2;
-   float  *fl,*fp,*mu,**jac,*ym1;
-   float  *fluxtop,fluxCH4;
-   float  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
-   float  rr,np,nm,factdec,s,test,time,ts,v;
-   double *fd,**jacd;
-   char   str2[15];
-   FILE   *out;
-
-/* va avec htoh2 */
-   float  dyh,dyh2;
-
-/* va avec aer */
-   float  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
-   float  **k_dep,**faer;
-   float  *productaer,*csurn,*csurh,*mmolaer;
-
-   if( (*aerprod) == 1 )
-   {
-    k_dep = rm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
-    faer  = rm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
-    productaer  = rm1d( 0, 3 );     /* local production rate by pathways */
-    mmolaer     = rm1d( 0, 3 );     /* local molar mass by pathways */
-    csurn = rm1d( 0, 3 );           /* local C/N by pathways */
-    csurh = rm1d( 0, 3 );           /* local C/H by pathways */
-   }
-
-/* DEBUG
-      printf("CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
-*/
-
-   for( i = 0; i <= NC; i++)
-   {
-     strcpy( corps[i], CORPS[i] );
-     corps[i][strcspn(CORPS[i], " ")] = '\0';
-   }
-   
-   strcpy( outlog, "chimietitan" );
-   strcat( outlog, ".log" );
-   deltamax = 1.e5;
-
-/*  DEBUG
-            out = fopen( outlog, "a" );
-   fprintf(out,"CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
-            fclose( out );
-*/
-   ym1     = rm1d( 0,   NC-1 );
-   fl      = rm1d( 0,   NC-1 );
-   fp      = rm1d( 0,   NC-1 );
-   fd      = dm1d( 0,   NC-1 );
-   mu      = rm1d( 0, NLEV-1 );
-   fluxtop = rm1d( 0,   NC-1 );
-   jac     = rm2d( 0,   NC-1, 0,   NC-1 );
-   jacd    = dm2d( 0,   NC-1, 0,   NC-1 );
-
-/* DEBUG */
-/*
-            out = fopen( "err.log", "a" );
-            fprintf( out,"%s\n", );
-            fclose( out );
-*/
-
-/* initialisation krate pour dissociations */
-
-   if( ( (*DECLIN) *10 + 267 ) < 14. ) 
-   {
-     declinint = 0;
-     dec = 0;
-   }
-   else 
-   {  
-       if( ( (*DECLIN) * 10 + 267 ) > 520. ) 
-       {
-          declinint = 14;
-          dec = 534;
-       }
-       else 
-       {
-          declinint = 1;
-          dec       = 27;
-          while( ( (*DECLIN) * 10 + 267 ) >= (float)(dec+20) ) 
-          {
-            dec += 40;
-            declinint++;
-          }
-       }
-   }
-   if( ( (*DECLIN) >= -24. ) && ( (*DECLIN) <= 24. ) ) 
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
-   else
-           factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;
-
-   for( i = 0; i <= RDISS; i++ )   /* RDISS correspond a la dissociation de N2 par les GCR */
-      for( j = 0; j <= NLEV-1; j++ )
-      { 
-         if( factdec < 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint-1] * fabs(factdec) 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 + factdec);
-         if( factdec > 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint+1] * factdec 
-                                         + KRPD[*LAT][j][i][declinint] * ( 1 - factdec);
-         if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
-      }
-        
-/* initialisation mu, CH4 au sol */
-     
-   for( j = 0; j <= NLEV-1; j++ )
-   {
-      mu[j] = 0.0e0;
-      for( i = 0; i <= ST-1; i++ ) 
-      {
-         if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
-         {
-             fluxCH4 = (*botCH4 - Y[i][j]);
-             Y[i][j] = *botCH4;
-         }
-         mu[j] += ( MASS[i] * Y[i][j] );
-      }
-   }
-
-/* ****************** */
-/*  Main loop: level  */
-/* ****************** */
-
-   for( j = NLEV-1; j >= 0; j-- )
-   {
-
-/* DEBUG
-            out = fopen( outlog, "a" );
-        fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
-            fclose( out );
-
-            out = fopen( "profils.log", "w" );
-    fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
-    for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
-    for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
-            fclose( out );
-*/
-
-      time     = ts = 0.0e0;
-/*      delta    = (*FIN);  */
-      delta    = 1.e-3;
-
-      for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
-
-/* ++++++++++++ */
-/*  time loop.  */
-/* ++++++++++++ */
-
-      while( time < (*FIN) )     
-      {
-
-/* Calcul de f et de la matrice jacobienne */
-/* --------------------------------------- */
-
-         for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
-         {
-            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
-
-            fp[i] = fl[i] = 0.0e0;                    /* init for next step */
-            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
-
-            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
-            {
-               ireac = prod[i][l];                  /* Number of the reaction involves. */
-               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
-               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
-               {
-                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
-                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
-               }
-               else                                 /* General case. */
-               {
-                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
-                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
-                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
-                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
-               }
-            }
-            
-            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
-            {
-               ireac = perte[i][l][0];              /* Reaction number. */
-               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
-               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
-               {
-                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
-                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
-               }
-               else                                 /* General case. */
-               {
-                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
-                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
-                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
-               }
-            }
-         }
-
-/* Aerosols */
-/* -------- */
-         if( (*aerprod) == 1 )
-         {
-             aer(corps,TEMP,NB,Y,&j,k_dep,faer,
-              &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
-              mmolaer,productaer,csurn,csurh);
-
-             for( i = 0; i <= 3; i++ )
-             {
-               PRODAER[i][j] = productaer[i];
-                  MAER[i][j] = mmolaer[i];
-                   CSN[i][j] = csurn[i];
-                   CSH[i][j] = csurh[i];
-             }
-/* DEBUG
-printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
-if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
-      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
-              fp[utilaer[2]],dyc2h2*NB[j]);
-if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
-      printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
-              fp[utilaer[5]],dyhcn*NB[j]);
-if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
-      printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
-              fp[utilaer[6]],dyhc3n*NB[j]);
-if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
-      printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
-              fp[utilaer[13]],dynccn*NB[j]);
-if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
-      printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
-              fp[utilaer[14]],dych3cn*NB[j]);
-if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
-      printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
-              fp[utilaer[15]],dyc2h3cn*NB[j]);
-*/
-
-             fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
-             fp[utilaer[5]] -= (    dyhcn * NB[j] );
-             fp[utilaer[6]] -= (   dyhc3n * NB[j] );
-             fp[utilaer[13]]-= (   dynccn * NB[j] );
-             fp[utilaer[14]]-= (  dych3cn * NB[j] );
-             fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
-             if( Y[utilaer[2]][j]  != 0.0 )
-       jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
-             if( Y[utilaer[5]][j]  != 0.0 )
-       jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
-             if( Y[utilaer[6]][j]  != 0.0 )
-       jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
-             if( Y[utilaer[13]][j] != 0.0 )
-     jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
-             if( Y[utilaer[14]][j] != 0.0 )
-     jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
-             if( Y[utilaer[15]][j] != 0.0 )
-     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
-         }
-            
-/* H -> H2 on haze particles */
-/* ------------------------- */
-         if( (*htoh2) == 1 )
-         {
-              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
-/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
-
-/* DEBUG
-printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
-if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
-printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
-if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
-printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
-*/
-
-              fp[utilaer[0]] += ( dyh  * NB[j] );
-              fp[utilaer[1]] += ( dyh2 * NB[j] );
-              if( Y[utilaer[0]][j] != 0.0 )
-       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
-         }
-
-         for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
-         {
-            for( k = 0; k <= ST-1; k++ )
-            {
-               jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
-               if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
-               jacd[i][k] = (double)jac[i][k];
-            }
-            
-            fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
-         }
-
-/*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
-         {
-            Y[i][j] = ??? ;
-         }
-*/
-
-/* Inversion of matrix cf method LU */
-/* -------------------------------- */
-
-/* inversion by blocs: */
-/* Hydrocarbons */
-
-         solve( jacd, fd, 0, NHC-1 );             
-
-         for( i = 0; i <= NHC-1; i++ )
-         {
-            Y[i][j] += (float)fd[i];
-            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
-         }
-
-/* Nitriles */
-
-         solve( jacd, fd, NHC, ST-1 );             
-
-         for( i = NHC+1; i <= ST-1; i++ )
-         {
-            Y[i][j] += (float)fd[i];
-            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
-         }
-
-/* end inversion by blocs: */
-
-         for( i = 0; i <= ST-1; i++ )
-         {
-
-/* CH4 au sol */
-/* ---------- */
-
-            if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
-            {
-                fluxCH4 += (*botCH4 - Y[i][0]);
-                Y[i][0] = *botCH4;
-            }
-
-         }
-         
-/* test evolution delta */
-/* -------------------- */
-
-         for( i = 0; i <= ST-1; i++ )
-         {
-            test = 1.0e-15; 
-            if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
-            {
-               conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
-
-               if( conv > ts )
-               {
-/*
-                  if( conv >= 0.1 )
-                  {
-                     out = fopen( outlog, "a" );
-                     fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
-                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
-                     fclose( out );
-                  }
-*/
-                  ts = conv;
-               }
-            }
-         }
-
-         if( ts < 0.1e0 )
-         {
-            for( i = 0; i <= ST-1; i++ )
-                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
-                 {
-                  out = fopen( outlog, "a" );
-                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
-                           corps[i], ym1[i], Y[i][j], j );
-                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
-                  fclose( out );
-                /*  exit(0); */
-                  Y[i][j] = 1.e-20;
-                 }
-            for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
-            time += delta;
-            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
-            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
-            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
-            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
-            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
-         
-            delta = min( deltamax, delta );
-         }
-         else
-         {
-            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
-
-            if(   ts > 0.8 )                    delta *= 1.e-6;
-            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
-            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
-            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
-            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
-            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
-         }
-         ts = 0.0e0;
-      }               
-
-/* +++++++++++++++++++ */
-/*  end of time loop.  */
-/* +++++++++++++++++++ */
-
-      for( i = 0; i <= ST-1; i++ ) 
-         if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
-            fluxCH4 *= ( MASS[i]/(6.022e23*time) ); 
-
-   }
-   
-/* **************** */        
-/* end of main loop */
-/* **************** */        
-
-/* Plafond: maintient de la composition */
-/* ------------------------------------ */
-/* 
-   for( i = 0; i <= ST-1; i++ ) 
-      if( YTOP[i] != 0.0e0 )
-      {
-             fluxtop[i]  = (Y[i][NLEV-1] - YTOP[i]) * MASS[i]/(6.022e23*(*FIN));
-             Y[i][NLEV-1]= YTOP[i];
-      }
-*/
-/* Plafond: !! OU !!  flux vertical     */
-/* ------------------------------------ */
-
-   for( i = 0; i <= ST-1; i++ ) 
-      if( FTOP[i] != 0.0e0 )
-      {
-             fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
-             Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
-             Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
-// on ajuste aussi le niveau dans la derniere couche...
-// pour eviter les effets vers le haut
-             Y[i][NLEV-1]  = Y[i][NLEV-2];
-      }
-
-/* Niveau de N2 */
-/* ------------ */
-      
-   for( j = 0; j <= NLEV-1; j++ ) 
-   {
-      conv = 0.0e0;
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
-      for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
-   }
-
-   if( (*aerprod) == 1 )
-   {
-    frm2d( k_dep, 1, 5, 1 );
-    frm2d(  faer, 1, 5, 1 );
-    frm1d( productaer,  0 );
-    frm1d( mmolaer,  0 );
-    frm1d( csurn, 0 );
-    frm1d( csurh, 0 );
-   }
-
-   frm1d(     ym1, 0 );
-   frm1d(      fl, 0 );
-   frm1d(      fp, 0 );
-   fdm1d(      fd, 0 );
-   frm1d(      mu, 0 );
-   frm1d( fluxtop, 0 );
-   frm2d(     jac, 0,   NC-1, 0 );
-   fdm2d(    jacd, 0,   NC-1, 0 );
-}
Index: trunk/LMDZ.TITAN/libf/chimtitan/omega.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/omega.c	(revision 1126)
+++ trunk/LMDZ.TITAN/libf/chimtitan/omega.c	(revision 1126)
@@ -0,0 +1,14 @@
+/* omega: correction to perfect gas compound */
+/* 30 Oct 96 */
+
+#include "titan.h"
+
+double omega( ts, epsa, epsb )
+double epsa,epsb,ts;
+{
+   double t;
+
+   t = ts / sqrt( epsa * epsb );
+   return 1.06036e0 * pow( t, -0.1561e0 ) + 0.193e0 * exp( -0.47635e0 * t )
+          + 1.03587e0 * exp( -1.52996e0 * t ) + 1.76474e0 * exp( -3.89411e0 * t );
+}
Index: trunk/LMDZ.TITAN/libf/chimtitan/solve.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/solve.c	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/chimtitan/solve.c	(revision 1126)
@@ -2,20 +2,22 @@
 /* cf Numerical Recipes LU Method for equations numbers */
 /* GCCM */
-/* une seule cellule concernee */
+/* similaire a inv (GP) */
 /* la matrice a est inversee seulement sur le bloc [n0;n1][n0;n1] */
-/* la matrice f ressort les dy */
 
 #include "titan.h"
 
-void solve( double **a, double *f, int n0, int n1 )
+void solve( double ***aa, int m, int n0, int n1 )
 {
-   int    i,ii,imax,k,l,ll;
-   double aamax,dum,*indx,sum,*vv,tmp;
+   int    i,ii,imax,j,k,l,ll;
+   double **a,aamax,dum,*indx,sum,**vv,tmp;
    FILE   *out;
 
    indx = dm1d( n0, n1 );
-   vv   = dm1d( n0, n1 );
+   vv   = dm2d( n0, n1, n0, n1 );
+   a    = dm2d( n0, n1, n0, n1 );
    imax = n0;
 
+   for( i = n0; i <= n1; i++ ) for( j = n0; j <= n1; j++ ) a[i][j]=aa[m][i][j];
+   
    for( i = n0; i <= n1; i++ )
    {
@@ -32,5 +34,5 @@
 /*         aamax = 1.e-30; */
       }
-      vv[i] = 1.0e0 / aamax;   /* Save the scaling */
+      vv[i][1] = 1.0e0 / aamax;   /* Save the scaling */
 /*
       if( (aamax > 1.0e100)||(aamax < 1.0e-100) )
@@ -59,5 +61,5 @@
             sum -= ( a[i][l] * a[l][k] );
          a[i][k] = sum;
-         dum        = vv[i] * fabs(sum);   /* Figure of merit for the pivot */
+         dum        = vv[i][1] * fabs(sum); /* Figure of merit for the pivot */
          if( dum >= aamax )          /* Is it better than the best so far ? */
          {
@@ -74,5 +76,5 @@
             a[k][l]    = dum;
          }
-         vv[imax] = vv[k];    /* Also interchange the scale factor */
+         vv[imax][1] = vv[k][1];    /* Also interchange the scale factor */
       }
       indx[k] = imax;
@@ -94,26 +96,41 @@
    }                                  /* Go back to the next column in the reduction */
 
+   for( i = n0; i <= n1; i++ )
+   {
+      for( k = n0; k <= n1; k++ )
+         vv[i][k] = 0.0e0;
+      vv[i][i] = 1.0e0;
+   }
+
+   for( l = n0; l <= n1; l++ )
+   {
       ii = n0-1;
       for( i = n0; i <= n1; i++ )
       {
          ll    = indx[i];
-         sum   = f[ll];
-         f[ll] = f[i];
+         sum   = vv[ll][l];
+         vv[ll][l] = vv[i][l];
          if( ii != (n0-1) )
             for( k = ii; k < i; k++ )
-               sum -= ( a[i][k] * f[k] );
+               sum -= ( a[i][k] * vv[k][l] );
          else if( sum != 0.0e0 ) ii = i;
-         f[i] = sum;
+         vv[i][l] = sum;
       }
       for( i = n1; i >= n0; i-- )
       {
-         sum = f[i];
+         sum = vv[i][l];
          if( i < n1 )
             for( k = i+1; k <= n1; k++ )
-               sum -= ( a[i][k] * f[k] );
-         f[i] = sum / a[i][i];
+               sum -= ( a[i][k] * vv[k][l] );
+         vv[i][l] = sum / a[i][i];
       }
+   }
          
+   for( i = n0; i <= n1; i++ ) 
+     for( l = n0; l <= n1; l++ ) 
+       aa[m][i][l]=vv[i][l];
+
    fdm1d( indx, n0 );
-   fdm1d(   vv, n0 );
+   fdm2d(   vv, n0, n1, n0 );
+   fdm2d(    a, n0, n1, n0 );
 }
Index: trunk/LMDZ.TITAN/libf/chimtitan/solve_b.c
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/solve_b.c	(revision 1126)
+++ trunk/LMDZ.TITAN/libf/chimtitan/solve_b.c	(revision 1126)
@@ -0,0 +1,126 @@
+/* matrix inversion */
+/* cf Numerical Recipes LU Method for equations numbers */
+/* GCCM */
+/* similaire a inv (GP) */
+/* la matrice a est inversee seulement sur le bloc [n0;n1][n0;n1] */
+
+#include "titan.h"
+
+void solve_b( double ***aa, double **f, int m, int n0, int n1 )
+{
+   int    i,ii,imax,j,k,l,ll;
+   double **a,aamax,dum,*indx,sum,*vv,tmp;
+   FILE   *out;
+
+   indx = dm1d( n0, n1 );
+   vv   = dm1d( n0, n1 );
+   a    = dm2d( n0, n1, n0, n1 );
+   imax = n0;
+
+   for( i = n0; i <= n1; i++ ) for( j = n0; j <= n1; j++ ) a[i][j]=aa[m][i][j];
+   
+   for( i = n0; i <= n1; i++ )
+   {
+      aamax = 0.0e0;
+      for( k = n0; k <= n1; k++ )
+         if( (tmp=fabs(a[i][k])) > aamax ) aamax = tmp;
+      if( aamax < 1.0e-20 )
+      {
+         out = fopen( "err.log", "a" );
+         fprintf( out, "Singular matrix. n0=%ld k=%ld aamax=%le\n",
+                        n0,k,aamax);
+         fclose( out );
+         exit(0);   
+/*         aamax = 1.e-30; */
+      }
+      vv[i] = 1.0e0 / aamax;   /* Save the scaling */
+/*
+      if( (aamax > 1.0e100)||(aamax < 1.0e-100) )
+      {
+         out = fopen( "err.log", "a" );
+         fprintf( out, "ATTENTION aamax = %le\n", aamax );
+         fclose( out );
+         exit( 0 );
+      }
+*/
+   }
+   for( k = n0; k <= n1; k++ )
+   {
+      for( i = n0; i < k; i++ )      /* This is equation 2.3.12 except for i = j */
+      {
+         sum = a[i][k];
+         for( l = n0; l < i; l++ )
+            sum -= ( a[i][l] * a[l][k] );
+         a[i][k] = sum;
+      }
+      aamax = 0.0e0;                 /* Initialize for the search for largest pivot element */
+      for( i = k; i <= n1; i++ )     /* This is i = j of equation 2.3.12 and */
+      {
+         sum = a[i][k];              /* i = J + 1,...,N of equation 2.3.13 */
+         for( l = n0; l < k; l++ )
+            sum -= ( a[i][l] * a[l][k] );
+         a[i][k] = sum;
+         dum        = vv[i] * fabs(sum); /* Figure of merit for the pivot */
+         if( dum >= aamax )          /* Is it better than the best so far ? */
+         {
+            imax  = i;
+            aamax = dum;
+         }
+      }
+      if( k != imax )                /* Do we need to interchange rows ? */
+      {
+         for( l = n0; l <= n1; l++ ) /* Yes, do so... */
+         {
+            dum        = a[imax][l];
+            a[imax][l] = a[k][l];
+            a[k][l]    = dum;
+         }
+         vv[imax] = vv[k];    /* Also interchange the scale factor */
+      }
+      indx[k] = imax;
+      if( fabs(a[k][k]) < 1.0e-20 )
+      {
+         out = fopen( "err.log", "a" );
+         fprintf( out, "Pivot too small. n0=%ld k=%ld fabs(a[k][k])=%le\n",
+                        n0,k,fabs(a[k][k]) );
+         fclose( out );
+         exit(0);   
+/*         a[k][k] = 1.e-20; */
+      }
+      if( k != n1 )                   /* If the pivot element is less than 1.0d-20 we */
+      {                               /* assume that the matrix is singular */
+          dum = a[k][k];              /* ( at least to the precision of the algorithm and the machine ) */
+         for( i = k+1; i <= n1; i++ ) /* Now, finally devide by the pivot element */
+            a[i][k] /= dum;
+      }
+   }                                  /* Go back to the next column in the reduction */
+
+   for( i = n0; i <= n1; i++ ) vv[i]=f[i][m];
+
+      ii = n0-1;
+      for( i = n0; i <= n1; i++ )
+      {
+         ll    = indx[i];
+         sum   = vv[ll];
+         vv[ll] = vv[i];
+         if( ii != (n0-1) )
+            for( k = ii; k < i; k++ )
+               sum -= ( a[i][k] * vv[k] );
+         else if( sum != 0.0e0 ) ii = i;
+         vv[i] = sum;
+      }
+      for( i = n1; i >= n0; i-- )
+      {
+         sum = vv[i];
+         if( i < n1 )
+            for( k = i+1; k <= n1; k++ )
+               sum -= ( a[i][k] * vv[k] );
+         vv[i] = sum / a[i][i];
+      }
+         
+   for( i = n0; i <= n1; i++ ) f[i][m]=vv[i];
+
+   fdm1d( indx, n0 );
+   fdm1d(   vv, n0 );
+   fdm2d(    a, n0, n1, n0 );
+}
Index: trunk/LMDZ.TITAN/libf/chimtitan/titan.h
===================================================================
--- trunk/LMDZ.TITAN/libf/chimtitan/titan.h	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/chimtitan/titan.h	(revision 1126)
@@ -7,11 +7,18 @@
 
 #define R0    (double)(2575.0) /* Titan's radius */
-#define NLEV  (int)(55)   /* Nbre de niv verticaux - aussi dans titan_for.h */
+#define NLEV  (int)(125)  /* Nbre de niv verticaux - =llm+70 dans common_mod */
+#define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff */
+#define NLRT  (int)(650)  /* Nbre de niv verticaux dans table fmoy - aussi dans common_mod */
+
+/* fluxes at 1300 km : upward is +, downward is - */
+#define top_H   (double)(+1.1e4)
+#define top_H2  (double)(+3.7e3)
+#define top_N4S (double)(-1.1e8) /* = -2.5e8/2.27 ... */
 
 /* DEPEND DE LA VERSION CHIMIE: */
 #define VERCHIM "chimie_simpnit_051006_bis"
-#define NREAC (int)(377)    /* nombre de reactions - aussi dans titan_for.h */
-#define RDISS (int)(54)     /* nombre de photodiss - aussi dans titan_for.h */
-#define NC    (int)(44)     /* nb de composes      - aussi dans titan_for.h */
+#define NREAC (int)(377)    /* nombre de reactions - aussi dans common_mod */
+#define RDISS (int)(54)     /* nombre de photodiss - aussi dans common_mod */
+#define NC    (int)(44)     /* nb de composes      - aussi dans common_mod */
 #define ST    (int)(NC)     /* nb de composes inverses */
 #define NHC   (int)(32)     /* nb hydrocarbons */
@@ -27,14 +34,11 @@
 #endif
 
-/* void  gptitan_(const int *, float *, float *, float *, float *, float *, 
-      char (*)[10], float (*)[NLEV], float *,
-      float *, float *, float *, int *, float *, 
-      float *, float (*)[NLEV][RDISS+1][15], float (*)[NLEV],
-      int (*)[5], int *, int *, int (*)[200], int (*)[200][2] ); */
 void  chimie_(char (*)[10], double *, double *, double (*)[NLEV], 
               int (*)[5], int *, int *, int (*)[200][2], int (*)[200]);
-void  comp_(char (*)[10], double *);
-void  disso_(double (*)[NLEV][RDISS+1][15], int *);
-void  solve( double **, double *, int, int );
+void  comp_(char (*)[10], double *, double *, double *, double (*)[NLEV]);
+void  disso_(double (*)[NLRT][RDISS+1][15], int *);
+double omega( double, double, double );
+void  solve( double ***, int, int, int );
+void  solve_b( double ***, double **, int, int, int );
 float *rm1d( int, int );
 float **rm2d( int, int, int, int );
Index: trunk/LMDZ.TITAN/libf/phytitan/calchim.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calchim.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/calchim.F	(revision 1126)
@@ -1,4 +1,4 @@
       SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
-     .                   ctemp,cplay,cplev,
+     .                   ctemp,cplay,cplev,czlay,czlev,
      .                   dqyc)
       
@@ -10,12 +10,13 @@
 c            adaptation pour Titan 3D: 02/2009
 c            adaptation pour // : 04/2013
-c
+c            extension chimie jusqu a 1300km : 10/2013
+c 
 c-------------------------------------------------
 c
       use dimphy
       use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
-     .                     NLEV,NC,ND,NR
+     .                     NLEV,NLRT,NC,ND,NR
       USE comgeomphy,  only: rlatd
-      use moyzon_mod,  only: klat
+      USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
       implicit none
 #include "dimensions.h"
@@ -35,5 +36,7 @@
       REAL         ctemp(nlon,klev)      ! Temperature
       REAL         cplay(nlon,klev)      ! pression (Pa)
-      REAL         cplev(nlon,klev)      ! pression intercouches (Pa)
+      REAL         cplev(nlon,klev+1)    ! pression intercouches (Pa)
+      REAL         czlay(nlon,klev)      ! altitude (m)
+      REAL         czlev(nlon,klev+1)    ! altitude intercouches (m)
       
       REAL         dqyc(nlon,klev,NC)    ! Tendances especes chimiques 
@@ -46,17 +49,19 @@
 c variables envoyees dans la chimie: double precision
 
-      REAL  temp_c(klev),press_c(klev)     ! T,p(mbar) a 1 lat donnee
-      REAL  declin_c                       ! declinaison en degres
-      REAL  cqy(klev,NC)                   ! frac mol qui seront modifiees
-      REAL  surfhaze(klev)
-      REAL  cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4)
-      REAL  rinter(klev),nb(klev)
+      REAL  temp_c(NLEV)
+      REAL  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
+      REAL  cqy(NLEV,NC)    ! frac mol qui seront modifiees
+      REAL  cqy0(NLEV,NC)    ! frac mol avant modif
+      REAL  surfhaze(NLEV)
+      REAL  cprodaer(NLEV,4),cmaer(NLEV,4)
+      REAL  ccsn(NLEV,4),ccsh(NLEV,4)
+! rmil: milieu de couche, grille pour K,D,p,ct,T,y
+! rinter: intercouche (grille RA dans la chimie)
+      REAL  rmil(NLEV),rinter(NLEV),nb(NLEV)
+      REAL,save :: kedd(NLEV)
+
       REAL  a,b
       character str1*1,str2*2
-      integer ntotftop
-      parameter (ntotftop=50)
-      integer nftop(ntotftop),isaisonflux
-      REAL  fluxtop(NC),latit,tabletmp(ntotftop),factflux
-      character*10 nomsftop(ntotftop+1)
+      REAL  latit
       character*20 formt1,formt2,formt3
       
@@ -67,14 +72,20 @@
       SAVE    firstcal
       
-      REAL  mass(NC),duree
-      REAL  tablefluxtop(NC,jjp1,5)
-      REAL  botCH4
+      integer dec,declinint,ialt
+      REAL  declin_c                       ! declinaison en degres
+      real  factalt,factdec,krpddec,krpddecp1,krpddecm1
+      real  duree
+      REAL,save :: mass(NC)
+      REAL,save,allocatable :: md(:,:)
+      REAL,save :: botCH4
       DATA  botCH4/0.05/
+      real,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
       REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
-      integer reactif(5,NR),nom_prod(NC),nom_perte(NC)
-      integer prod(200,NC),perte(2,200,NC)
-      SAVE    mass,tablefluxtop,botCH4
-      SAVE    reactif,nom_prod,nom_perte,prod,perte
-      
+      integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
+      integer,save :: prod(200,NC),perte(2,200,NC)
+      character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
+      real  ficalt,oldq,newq,zalt
+      logical okfic
+
 c-----------------------------------------------------------------------
 c***********************************************************************
@@ -94,16 +105,7 @@
 c ************************************
 
-        allocate(krpd(15,ND+1,klev,jjp1),krate(klev,NR))
-
-c Verification dimension verticale: coherence titan_for.h et klev
-c --------------------------------
-
-        if (klev.ne.NLEV) then
-	   print*,'PROBLEME de coherence dans la dimension verticale:'
-     .           ,klev,NLEV
-           STOP
-	endif
-
-c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro
+        allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
+
+c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
 c ----------------------------------
 
@@ -114,114 +116,72 @@
 	endif
 
-c calcul de temp_c, densites et press_c au milieu de l'ensemble des points:
-c ----------------------------------------------------------------------
-
-        print*,'pression, densites et temp (chimie):'
+c calcul de temp_c, densites et press_c en moyenne planetaire :
+c --------------------------------------------------------------
+
+        print*,'pression, densites et temp (init chimie):'
         print*,'level, press_c, nb, temp_c'
         DO l=1,klev
+         rinter(l)  = (zlevmoy(l)+RA)/1000.
+         rmil(l)    = (zlaymoy(l)+RA)/1000.
 c     temp_c (K):
-         temp_c(l)  = ctemp(nlon/2+1,l)
+         temp_c(l)  = tmoy(l)
 c     press_c (mbar):
-         press_c(l) = cplay(nlon/2+1,l)/100.
+         press_c(l) = playmoy(l)/100.
 c     nb (cm-3):
          nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
-         print*,l,press_c(l),nb(l),temp_c(l)
+         print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
         ENDDO
+
+c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+1,NLEV
+         rinter(l) = rinter(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+         rmil(l)   = rmil(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+       enddo        
+
+c lecture de tcp.ver, une seule fois
+c remplissage r1d,t1d,ct1d,p1d
+        open(11,file='../../INPUT/tcp.ver',status='old')
+        read(11,*) dummy
+        do i=1,131
+          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
+c         print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
+        enddo 
+        close(11)
+
+c extension pour klev+1 a NLEV avec tcp.ver
+c la jonction klev/klev+1 est brutale... Tant pis ?
+        ialt=1
+        do l=klev+1,NLEV
+           zalt=rmil(l)-RA/1000.
+           do i=ialt,130
+            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
+              ialt=i
+            endif
+           enddo
+           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
+           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt) 
+     &                      + log(p1d(ialt+1)) * factalt    )
+           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) 
+     &                      + log(ct1d(ialt+1)) * factalt    )
+           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
+           print*,l,zalt,press_c(l),nb(l),temp_c(l)
+        enddo
         
 c caracteristiques des composes:        
 c -----------------------------
-        mass = 0.0
-        call comp(nomqy_c,mass)
+        mass(:) = 0.0
+        call comp(nomqy_c,nb,temp_c,mass,md)
         print*,'           Mass'
         do ic=1,NC
           print*,nomqy_c(ic),mass(ic)
+c         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
         enddo
 	
 
-c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON... 
-
-c flux dans la derniere couche:(issu du modele 1D, a 500 km) 
-c -----------------------------
-c       !!  en cm-2.s-1 !!
-c  ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54 
-c 
-c Lecture des tables de flux en fonction lat et saison
-c ----------------------------------------------------
-
-	WRITE(str2(1:2),'(i2.2)') ntotftop
-        formt1 = '(A10,'//str2//'(3X,A10))'
-        formt2 = '(F12.6,'//str2//'(2X,F13.6))'
-        formt3 = '(F13.6,'//str2//'(2X,F13.6))'
-	
-	do j=1,jjp1
-	  do ic=1,NC
-	    do l=1,5
-	     tablefluxtop(ic,j,l) = 0.
-	    enddo
-	  enddo
-	enddo
-	
-        do l=1,4
-	  WRITE(str1(1:1),'(i1.1)') l
-c hx1 -> Ls=RPI/2 / hx4 -> Ls=0
-          open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1)
-c         open(10,file="FLUXTOP/flux500.hs"//str1)
-c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ...
-
-          read(10,formt1) nomsftop
-	  do j=1,ntotftop
-           do i=1,10   !justification a gauche...
-             if (nomsftop(j+1)(i:i).ne.' ') then
-                nomsftop(j) = nomsftop(j+1)(i:)
-                goto 100
-             endif
-           enddo
-100        continue
-c         print*,nomsftop(j)
-           nftop(j) = 0
-           do i=1,NC
-	    if (nomqy_c(i).eq.nomsftop(j)) then
-             nftop(j) = i
-	    endif
-	   enddo
-c          if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j))
-	  enddo
-
-	 
-c lecture des flux. Les formats sont un peu alambiques... 
-c     a revoir eventuellement
-	  do j=1,jjp1/2+1
-	    read(10,formt2) latit,tabletmp
-	    do i=1,ntotftop
-              if (nftop(i).ne.0) then
-	       tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
-              endif
-	    enddo
-	  enddo
-	  do j=jjp1/2+2,jjp1
-	    read(10,formt3) latit,tabletmp
-	    do i=1,ntotftop
-              if (nftop(i).ne.0) then
-	       tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
-              endif
-	    enddo
-	  enddo
-
-          close(10)
-	enddo  ! l
-	
-	do j=1,jjp1
-	  do ic=1,NC
-	     tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5)
-	  enddo
-	enddo
-	
-c  Stockage des composes utilises dans la prod d'aerosols
+c  Stockage des composes utilises dans la prod d aerosols
 c     (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
 c     ! decalage de 1 car utilise dans le c !
-c ------------------------------------------
-c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN
-c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10
-c ------------------------------------------
 
         do ic=1,NC
@@ -244,9 +204,4 @@
 	  if (nomqy_c(ic).eq."HCN") then
             utilaer(6)  = ic-1
-	    do j=1,jjp1
-              do i=1,5
-             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6.
-              enddo
-            enddo
 	  endif
 	  if (nomqy_c(ic).eq."HC3N") then
@@ -268,10 +223,4 @@
 	  if (nomqy_c(ic).eq."C2H2") then
             utilaer(3)  = ic-1
-	    do j=1,jjp1
-              do i=1,5
-c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3
-             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10.
-              enddo
-            enddo
 	  endif
 	  if (nomqy_c(ic).eq."AC6H6") then
@@ -298,28 +247,4 @@
 	  if (nomqy_c(ic).eq."AC6H5") then
             utilaer(12) = ic-1
-	  endif
-
-	  if (nomqy_c(ic).eq."C2H6") then
-	    do j=1,jjp1
-              do i=1,5
-c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640.
-             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200.
-              enddo
-            enddo
-	  endif
-	  if (nomqy_c(ic).eq."C3H8") then
-	    do j=1,jjp1
-              do i=1,5
-c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
-             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.)
-              enddo
-            enddo
-	  endif
-	  if (nomqy_c(ic).eq."C4H10") then
-	    do j=1,jjp1
-              do i=1,5
-             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
-              enddo
-            enddo
 	  endif
 
@@ -369,4 +294,29 @@
 c       enddo
 
+
+c   init kedd
+c   ---------
+c   kedd stays constant with time and space 
+c   => linked to pressure rather than altitude...
+ 
+      kedd(:) = 5e2 ! valeur mise par defaut  
+               ! UNITE ? doit etre ok pour gptitan
+      do l=1,NLEV
+       zalt=rmil(l)-RA/1000.  ! z en km
+       if     ((zalt.ge.250.).and.(zalt.le.400.)) then 
+         kedd(l) = 10.**(3.+(zalt-250.)/50.)
+! 1E3 at 250 km
+! 1E6 at 400 km
+       elseif ((zalt.gt.400.).and.(zalt.le.600.)) then 
+         kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.)
+! 2E7 at 600 km
+       elseif ((zalt.gt.600.).and.(zalt.le.900.)) then 
+         kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.)
+! 1E8 above 900 km
+       elseif ( zalt.gt.900.                    ) then 
+        kedd(l) = 1.e8
+       endif
+      enddo
+
       ENDIF  ! premier appel
 
@@ -380,16 +330,4 @@
 c      print*,'declinaison en degre=',declin_c
        
-c-----------------------------------------------------------------------
-c
-c   calcul du facteur d'interpolation entre deux saisons pour flux
-c   --------------------------------------------------------------
-
-       isaisonflux = int(ls_rad*2./RPI)+1
-       if (isaisonflux.eq.5) then
-         isaisonflux = 1
-       endif
-       factflux =    (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/
-     .   (sin(isaisonflux*RPI/2.)-sin((isaisonflux-1)*RPI/2.))
-
 c***********************************************************************
 c***********************************************************************
@@ -410,4 +348,23 @@
 c***********************************************************************
 
+c-----------------------------------------------------------------------
+c
+c   Distance radiale (intercouches, en km)
+c   ----------------------------------------
+
+       do l=1,klev
+         rinter(l) = (RA+czlev(j,l))/1000.
+         rmil(l)   = (RA+czlay(j,l))/1000.
+c        print*,rinter(l)
+       enddo
+
+c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+1,NLEV
+         rinter(l) = rinter(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+         rmil(l)   = rmil(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+       enddo
+        
 c-----------------------------------------------------------------------
 c
@@ -423,51 +380,75 @@
                nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
        ENDDO
-       
-c-----------------------------------------------------------------------
-c
-c   Distances radiales (intercouches, en km)
-c   ----------------------------------------
-
-       rinter(1) = RA/1000.  
-       do l=2,klev
-c A REVOIR: g doit varier avec r !
-         rinter(l) = rinter(l-1) + 
-     .    (cplev(j,l-1)-cplev(j,l))/cplay(j,l)*(RD*temp_c(l)/RG)/1000.
-c        print*,rinter(l)
-       enddo
-
-c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-c FLUX AU TOP: interpolation en saison et
-c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C) 
-c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        do ic=1,NC
-          fluxtop(ic) = tablefluxtop(ic,j,isaisonflux)  *(1-factflux)
-     .                + tablefluxtop(ic,j,isaisonflux+1)*  factflux
-          fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5)
-        enddo
-
-c TEST: sortie de l'un des flux avec saveLs et factflux
-c       dans un fichier special pour tracer avec gnuplot
-c       if (j.eq.2) then
-c       open(unit=11,file="flux_80N.txt",status='old',position='append')
-c         write(11,*) ls_rad*180./RPI, factflux, 
-c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
-c         write(11,*) ls_rad*180./RPI, factflux, 
-c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
-c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
-c       close(11)
-c       endif
-c       if (j.eq.17) then
-c       open(unit=11,file="flux_eq.txt",status='old',position='append')
-c         write(11,*) ls_rad*180./RPI, factflux, 
-c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
-c         write(11,*) ls_rad*180./RPI, factflux, 
-c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
-c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
-c       close(11)
-c       endif
-c FIN TEST: sortie de l'un des flux avec ls et factflux
-
-c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+c extension pour klev+1 a NLEV avec tcp.ver
+       ialt=1
+       do l=klev+1,NLEV
+           zalt=rmil(l)-RA/1000.
+           do i=ialt,130
+            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
+              ialt=i
+            endif
+           enddo
+           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
+           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt) 
+     &                      + log(p1d(ialt+1)) * factalt    )
+           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) 
+     &                      + log(ct1d(ialt+1)) * factalt    )
+           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
+       enddo
+               
+c-----------------------------------------------------------------------
+c
+c   passage krpd => krate 
+c   ---------------------
+c initialisation krate pour dissociations 
+
+      if ((declin_c*10+267).lt.14.) then
+          declinint = 0
+          dec       = 0
+      else 
+       if ((declin_c*10+267).gt.520.) then
+          declinint = 14
+          dec       = 534
+       else 
+          declinint = 1
+          dec       = 27
+          do while( (declin_c*10+267).ge.real(dec+20) )
+            dec       = dec+40
+            declinint = declinint+1
+          enddo
+       endif
+      endif
+      if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
+          factdec = ( declin_c - (dec-267)/10. ) / 4.
+      else
+          factdec = ( declin_c - (dec-267)/10. ) / 2.7
+      endif
+
+      do l=1,NLEV
+
+c INTERPOL EN ALT POUR k (krpd tous les deux km)
+        ialt    = int((rmil(l)-RA/1000.)/2.)+1
+        factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
+
+        do i=1,ND+1
+          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
+          krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
+          krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
+
+                    ! ND+1 correspond a la dissociation de N2 par les GCR
+          if ( factdec.lt.0. ) then 
+        krate(l,i) = krpddecm1 * abs(factdec)
+     &             + krpddec   * ( 1 + factdec)
+          endif
+          if ( factdec.gt.0. ) then
+        krate(l,i) = krpddecp1 * factdec
+     &             + krpddec   * ( 1 - factdec)
+          endif
+          if ( factdec.eq.0. ) krate(l,i) = krpddec
+        enddo        
+      enddo        
 
 c-----------------------------------------------------------------------
@@ -475,11 +456,52 @@
 c   composition 
 c   ------------
-       
+
        do ic=1,NC
         do l=1,klev
-          cqy(l,ic) = qy_c(j,l,ic) 
-        enddo
-       enddo
-              
+          cqy(l,ic)  = qy_c(j,l,ic) 
+          cqy0(l,ic) = cqy(l,ic)
+        enddo
+       enddo
+
+c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
+
+      WRITE(str2,'(i2.2)') klat(j)
+      fich = "comp_"//str2//".dat"
+      inquire (file=fich,exist=okfic)
+      if (okfic) then
+       open(11,file=fich,status='old')
+c premiere ligne=declin
+       read(11,'(A15)') ficdec
+       write(curdec,'(E15.5)') declin_c
+c si la declin est la meme: ce fichier a deja ete reecrit 
+c on lit la colonne t-1 au lieu de la colonne t
+c (cas d une bande de latitude partagee par 2 procs)
+       do ic=1,NC
+        read(11,'(A10)') name
+        if (name.ne.nomqy_c(ic)) then
+          print*,"probleme lecture ",fich
+          print*,name,nomqy_c(ic)
+          stop
+        endif
+        if (ficdec.eq.curdec) then
+          do l=klev+1,NLEV
+            read(11,*) ficalt,cqy(l,ic),newq
+          enddo
+        else
+          do l=klev+1,NLEV
+            read(11,*) ficalt,oldq,cqy(l,ic)
+          enddo
+        endif
+       enddo
+       close(11)
+      else       ! le fichier n'est pas la
+       do ic=1,NC
+        do l=klev+1,NLEV
+          cqy(l,ic)=cqy(klev,ic)
+        enddo
+       enddo
+      endif
+      cqy0 = cqy
+
 c-----------------------------------------------------------------------
 c
@@ -502,7 +524,7 @@
        
        call gptitan(rinter,temp_c,nb,
-     $              nomqy_c,cqy,fluxtop,
-     $              declin_c,duree,(klat(j)-1),mass,
-     $              botCH4,krpd,krate,reactif,
+     $              nomqy_c,cqy,
+     $              duree,(klat(j)-1),mass,md,
+     $              kedd,botCH4,krate,reactif,
      $              nom_prod,nom_perte,prod,perte,
      $              aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
@@ -514,5 +536,5 @@
        do ic=1,NC
         do l=1,klev
-          dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim  ! en /s
+          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
         enddo
        enddo
@@ -535,4 +557,23 @@
 
        endif
+
+c-----------------------------------------------------------------------
+c
+c   sauvegarde compo sur NLEV 
+c   -------------------------
+
+c dans fichier compo_klat(j) (01 à 49)
+       
+      open(11,file=fich,status='replace')
+c premiere ligne=declin
+      write(11,'(E15.5)') declin_c
+      do ic=1,NC
+       write(11,'(A10)') nomqy_c(ic)
+       do l=klev+1,NLEV
+        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
+     .                                cqy0(l,ic),cqy(l,ic)
+       enddo
+      enddo
+      close(11)
 
 c***********************************************************************
Index: trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90	(revision 1126)
+++ trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90	(revision 1126)
@@ -0,0 +1,38 @@
+subroutine cirs_haze(press,wno,taeros,taeroscat)
+IMPLICIT NONE
+
+real,intent(in)   :: press,wno
+real,intent(inout):: taeros,taeroscat
+
+!---------------------------
+! Attention: taeros est en m-1
+! mais la valeur semble devoir etre multipliee par lambda (en m) pour
+! obtenir une valeur comparable a celle fournies par Sandrine...
+! a tirer au clair...
+!---------------------------
+
+real         :: taerosold
+logical,save :: firstcall=.true.
+
+if (firstcall) then
+   print*,"CIRS HAZE"
+   firstcall=.false.
+endif
+
+if (wno.eq.600.) then
+ print*,press,wno,taeros,taeroscat
+endif
+
+taerosold = taeros
+
+! modif de taeros
+
+! taeroscat est modifie proportionnellement a taeros
+
+if (taerosold.ne.0.) then
+   taeroscat = taeroscat/taerosold*taeros
+endif
+
+! Je maintiens le cbar du calcul microphysique
+
+end subroutine cirs_haze
Index: trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90	(revision 1126)
@@ -34,7 +34,9 @@
 
 ! ancien titan_for.h
-      INTEGER, PARAMETER :: NLEV=llm,NC=44,ND=54,NR=377
+      INTEGER, PARAMETER :: NLEV=llm+70,NC=44,ND=54,NR=377
 !$OMP THREADPRIVATE(NLEV,NC,ND,NR)
 !!!  doivent etre en accord avec titan.h
+! pour l'UV (650 niveaux de 2 km)
+      INTEGER, PARAMETER :: NLRT=650
 
 ! ancien diagmuphy.h
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F	(revision 1126)
@@ -162,4 +162,7 @@
           XICLDI2(K)=TNI
  420  CONTINUE
+
+c DEBUG
+c       print*,"wnoi=",WNOI
 
 C
@@ -284,4 +287,20 @@
         TauGID(ig,:,:) = TAUGID_1pt(:,:) 
 
+c DEBUG
+c     if(ig.eq.(ngrid/2+16)) then
+c         print*,ig,'/',KLON,':'
+c         print*,'TauHID_1',TAUHID(ig,1,:)
+c         print*,'TauGID_1',TAUGID(ig,1,:)
+c         print*,'TauHID_50',TAUHID(ig,50,:)
+c         print*,'TauGID_50',TAUGID(ig,50,:)
+c         print*,"DTAUI_1=",DTAUI(ig,1,:)
+c         print*,"DTAUI_50=",DTAUI(ig,50,:)
+c         print*,'cosBI_1',COSBI(ig,1,:)
+c         print*,'cosBI_50',COSBI(ig,50,:)
+c         print*,'WBARI_1',WBARI(ig,1,:)
+c         print*,'WBARI_50',WBARI(ig,50,:)
+c         stop
+c     endif
+
 c************************************************************************
 c************************************************************************
Index: trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F	(revision 1126)
@@ -282,4 +282,6 @@
 c     CBAR=xv3(j,k)
 
+c     print*, 'HERE, CIRS AEROSOLS'
+c     call cirs_haze(PRESS(J),WNOI(K),TAEROS,TAEROSCAT,CBAR)
 
 c*********** EN TRAVAUX ***************************
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F	(revision 1126)
@@ -148,4 +148,7 @@
 c       close(1)
 
+c DEBUG
+c       print*,"wnov=",WNOV
+
       endif    ! fin initialisations premier appel
 
@@ -167,5 +170,5 @@
         else
          if (ig.eq.1)  then
-c initialisation zqaer_1pt a partir d'une look-up table (uniforme en ig)
+c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
 c boucle sur nrad=10
            open(10,file="qaer_eq_1d.dat")
@@ -220,4 +223,14 @@
         TauGVD(ig,:,:) = TAUGVD_1pt(:,:) 
 
+c DEBUG
+c     if(ig.eq.(ngrid/2+16)) then
+c         print*,ig,'/',KLON,':'
+c         print*,'TauHVD_1',TAUHVD(ig,1,:)
+c         print*,'TauGVD_1',TAUGVD(ig,1,:)
+c         print*,'TauHVD_50',TAUHVD(ig,50,:)
+c         print*,'TauGVD_50',TAUGVD(ig,50,:)
+c     stop
+c     endif
+
  101  CONTINUE
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F	(revision 1126)
@@ -231,4 +231,7 @@
 
        IF (TAEROSCAT.le.0.) CBAR=0. 
+
+c     print*, 'HERE, CIRS AEROSOLS'
+c     call cirs_haze(PRESS(J),WNOV(K),TAEROS,TAEROSCAT,CBAR)
 
 1699  FORMAT(a3,2I3,3(ES15.7,1X))
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1126)
@@ -419,5 +419,5 @@
          itap    = 0
          itaprad = 0
-         itapchim    = 0
+         itapchim    = 1
 
 c init rnuabar
@@ -1126,5 +1126,9 @@
          tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
      .                        + d_tr_mph(:,:,1:nmicro)*dtime
-
+c        call WriteField_phy('physiq_d_tr_mph01',
+c    .                          d_tr_mph(:,:,1),klev)
+c        call WriteField_phy('physiq_d_tr_mph10',
+c    .                          d_tr_mph(:,:,10),klev)
+	endif
 c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
 c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
@@ -1141,6 +1145,4 @@
         ENDDO
 
-	endif
-
 c condensation:
 c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
@@ -1149,10 +1151,11 @@
      .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
         endif
+
+c chimie:
         if ((chimi).and.(nqmax.gt.nmicro)) then
-c chimie:
          tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
 	endif
 
-      endif
+      endif  !iflag_trac=1
 
 c       ch4=0.
Index: trunk/LMDZ.TITAN/libf/phytitan/phytrac.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phytrac.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/phytrac.F	(revision 1126)
@@ -459,5 +459,5 @@
 c ------------
        CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
-     .              ztemp,zplay,zplev,
+     .              ztemp,zplay,zplev,zzlay,zzlev,
      .              pdyfi)   
 c ychim ne doit pas etre modifie, pdyfi en /s
Index: trunk/LMDZ.TITAN/libf/phytitan/profile.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 1126)
@@ -94,7 +94,8 @@
        b2 =       3183. 
        c2 =       4737. 
+c pres est en Pa => conversion car expression veut p en hPa
        DO il=1,nlev
-         temp(il)=a1*exp(-((pres(il)-b1)/c1)**2.) 
-     .          + a2*exp(-((pres(il)-b2)/c2)**2.)
+         temp(il)=a1*exp(-((pres(il)/100.-b1)/c1)**2.) 
+     .          + a2*exp(-((pres(il)/100.-b2)/c2)**2.)
        ENDDO
        zkm(1)  = 0.0
