Ignore:
Timestamp:
Dec 17, 2013, 1:02:44 PM (11 years ago)
Author:
slebonnois
Message:

SL: update of Titan photochemical module to include computation of chemistry up to 1300 km

Location:
trunk/LMDZ.TITAN/libf/chimtitan
Files:
2 added
2 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/chimtitan/comp.c

    r3 r1126  
    44#include "titan.h"
    55
    6 void comp_(char CORPS[][10], double *MASS)
     6void comp_(char CORPS[][10], double *CT, double *TEMP,
     7           double *MASS, double MD[][NLEV])
    78{
    8    int   i;
    9    char  corps[100][10];
     9   int   i,j;
     10   char  corps[NC+1][10];
     11
     12   double m,ma,epsa,sig,siga,p;
     13
     14   p         = 2.976e07;   /*9 10^7 R / 8 pi */
     15
     16      /* WARNING BACKGROUND GAS IS N2 */
     17
     18   ma   = 28.0134e0;               /* mass of background gas in g */
     19   siga = 3.798e0;         /* Lennard-Jones length of background gas 1/10 nm */
     20   epsa = 71.4e0;       /* Lennard-Jones energy of background gas */
    1021
    1122   for( i = 0; i <= NC; i++)
    … …  
    1526   }
    1627
     28   for( i = 0; i <= NC-1; i++)
     29   {
     30      for( j = 0; j <= NLEV-1; j++ ) MD[i][j] = 0.0e0;
     31   }
     32
    1733   for( i = 0; i <= NC-1; i++ )
    1834   {
    … …  
    2036      {
    2137         MASS[i] = 16.04e0;
     38         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     39         sig     = 1.0e-16 * pow( ( siga + 3.758e0 ), 2.0e0 );
     40         for( j = 0; j <= NLEV-1; j++ )
     41            MD[i][j] = sqrt( p * TEMP[j] * m )
     42           / ( CT[j] * sig * omega(TEMP[j],epsa,148.6e0) );
    2243      }
    2344      if( strcmp(corps[i], "H") == 0 )
    … …  
    2849      {
    2950         MASS[i] = 2.0158e0;
     51         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     52         sig     = 1.0e-16 * pow( ( siga + 2.827e0 ), 2.0e0 );
     53         for( j = 0; j <= NLEV-1; j++ )
     54            MD[i][j] = sqrt( p * TEMP[j] * m )
     55           / ( CT[j] * sig * omega(TEMP[j],epsa,59.7e0) );
    3056      }
    3157      if( strcmp(corps[i], "CH") == 0 )
    3258      {
    3359         MASS[i] = 13.02e0;
     60         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     61         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
     62         for( j = 0; j <= NLEV-1; j++ )
     63            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    3464      }
    3565      if( ( strcmp( corps[i], "CH2" ) == 0 ) || ( strcmp( corps[i], "CH2s" ) == 0 ) )
    3666      {
    3767         MASS[i] = 14.03e0;
     68         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     69         sig     = 1.0e-16 * pow( ( siga + 3.4e0 ), 2.0e0 );
     70         for( j = 0; j <= NLEV-1; j++ )
     71            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    3872      }
    3973      if( strcmp(corps[i], "CH3") == 0 )
    4074      {
    4175         MASS[i] = 15.03e0;
     76         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     77         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
     78         for( j = 0; j <= NLEV-1; j++ )
     79            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    4280      }
    4381      if( strcmp(corps[i], "C") == 0 )
    4482      {
    4583         MASS[i] = 12.01e0;
     84         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     85         sig     = 1.0e-16 * pow( ( siga + 1.4e0 ), 2.0e0 );
     86         for( j = 0; j <= NLEV-1; j++ )
     87            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    4688      }
    4789      if( strcmp(corps[i], "C2") == 0 )
    4890      {
    4991         MASS[i] = 24.02e0;
     92         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     93         sig     = 1.0e-16 * pow( ( siga + 3.2e0 ), 2.0e0 );
     94         for( j = 0; j <= NLEV-1; j++ )
     95            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    5096      }
    5197      if( strcmp(corps[i], "C2H") == 0 )
    5298      {
    5399         MASS[i] = 25.03e0;
     100         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     101         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
     102         for( j = 0; j <= NLEV-1; j++ )
     103            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    54104      }
    55105      if( strcmp(corps[i], "C2H3") == 0 )
    56106      {
    57107         MASS[i] = 27.05e0;
     108         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     109         sig     = 1.0e-16 * pow( ( siga + 3.8e0 ), 2.0e0 );
     110         for( j = 0; j <= NLEV-1; j++ )
     111            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    58112      }
    59113      if( strcmp(corps[i], "C2H4") == 0 )
    60114      {
    61115         MASS[i] = 28.05e0;
     116         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     117         sig     = 1.0e-16 * pow( ( siga + 4.163e0 ), 2.0e0 );
     118         for( j = 0; j <= NLEV-1; j++ )
     119            MD[i][j] = sqrt( p * TEMP[j] * m )
     120           / ( CT[j] * sig * omega(TEMP[j],epsa,224.7e0) );
    62121      }
    63122      if( strcmp(corps[i], "C2H2") == 0 )
    64123      {
    65124         MASS[i] = 26.04e0;
     125         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     126         sig     = 1.0e-16 * pow( ( siga + 4.033e0 ), 2.0e0 );
     127         for( j = 0; j <= NLEV-1; j++ )
     128            MD[i][j] = sqrt( p * TEMP[j] * m )
     129           / ( CT[j] * sig * omega(TEMP[j],epsa,231.8e0) );
    66130      }
    67131      if( strcmp(corps[i], "C2H5") == 0 )
    68132      {
    69133         MASS[i] = 29.06e0;
     134         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     135         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     136         for( j = 0; j <= NLEV-1; j++ )
     137            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    70138      }
    71139      if( strcmp(corps[i], "C2H6") == 0 )
    72140      {
    73141         MASS[i] = 30.07e0;
     142         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     143         sig     = 1.0e-16 * pow( ( siga + 4.443e0 ), 2.0e0 );
     144         for( j = 0; j <= NLEV-1; j++ )
     145            MD[i][j] = sqrt( p * TEMP[j] * m )
     146           / ( CT[j] * sig * omega(TEMP[j],epsa,215.7e0) );
    74147      }
    75148      if( strcmp(corps[i], "C3H2") == 0 )
    76149      {
    77150         MASS[i] = 38.05e0;
     151         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     152         sig     = 1.0e-16 * pow( ( siga + 4.6e0 ), 2.0e0 );
     153         for( j = 0; j <= NLEV-1; j++ )
     154            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    78155      }
    79156      if( strcmp(corps[i], "C3H3") == 0 )
    80157      {
    81158         MASS[i] = 39.06e0;
     159         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     160         sig     = 1.0e-16 * pow( ( siga + 4.7e0 ), 2.0e0 );
     161         for( j = 0; j <= NLEV-1; j++ )
     162            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    82163      }
    83164      if( ( strcmp(corps[i], "CH2CCH2") == 0 ) || ( strcmp(corps[i], "CH3CCH") == 0 ) )
    84165      {
    85166         MASS[i] = 40.07e0;
     167         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     168         sig     = 1.0e-16 * pow( ( siga + 4.761e0 ), 2.0e0 );
     169         for( j = 0; j <= NLEV-1; j++ )
     170            MD[i][j] = sqrt( p * TEMP[j] * m )
     171           / ( CT[j] * sig * omega(TEMP[j],epsa,251.8e0) );
    86172      }
    87173      if( strcmp(corps[i], "C3H5") == 0 )
    88174      {
    89175         MASS[i] = 41.07e0;
     176         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     177         sig     = 1.0e-16 * pow( ( siga + 4.78e0 ), 2.0e0 );
     178         for( j = 0; j <= NLEV-1; j++ )
     179            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    90180      }
    91181      if( strcmp(corps[i], "C3H6") == 0 )
    92182      {
    93183         MASS[i] = 42.08e0;
     184         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     185         sig     = 1.0e-16 * pow( ( siga + 4.807e0 ), 2.0e0 );
     186         for( j = 0; j <= NLEV-1; j++ )
     187            MD[i][j] = sqrt( p * TEMP[j] * m )
     188           / ( CT[j] * sig * omega(TEMP[j],epsa,248.9e0) );
    94189      }
    95190      if( strcmp(corps[i], "C3H7") == 0 )
    96191      {
    97192         MASS[i] = 43.09e0;
     193         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     194         sig     = 1.0e-16 * pow( ( siga + 5.0e0 ), 2.0e0 );
     195         for( j = 0; j <= NLEV-1; j++ )
     196            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    98197       }
    99198      if( strcmp(corps[i], "C3H8") == 0 )
    100199      {
    101200         MASS[i] = 44.11e0;
     201         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     202         sig     = 1.0e-16 * pow( ( siga + 5.118e0 ), 2.0e0 );
     203         for( j = 0; j <= NLEV-1; j++ )
     204            MD[i][j] = sqrt( p * TEMP[j] * m )
     205           / ( CT[j] * sig * omega(TEMP[j],epsa,237.1e0) );
    102206      }
    103207      if( strcmp(corps[i], "C4H") == 0 )
    104208      {
    105209         MASS[i] = 49.05e0;
     210         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     211         sig     = 1.0e-16 * pow( ( siga + 4.2e0 ), 2.0e0 );
     212         for( j = 0; j <= NLEV-1; j++ )
     213            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    106214      }
    107215      if( ( strcmp(corps[i], "C4H2") == 0 )||( strcmp(corps[i], "C4H2s") == 0 ) )
    108216      {
    109217         MASS[i] = 50.06e0;
     218         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     219         sig     = 1.0e-16 * pow( ( siga + 4.3e0 ), 2.0e0 );
     220         for( j = 0; j <= NLEV-1; j++ )
     221            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    110222      }
    111223      if( strcmp(corps[i], "C4H3") == 0 )
    112224      {
    113225         MASS[i] = 51.07e0;
     226         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     227         sig     = 1.0e-16 * pow( ( siga + 4.4e0 ), 2.0e0 );
     228         for( j = 0; j <= NLEV-1; j++ )
     229            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    114230      }
    115231      if( strcmp(corps[i], "C4H4") == 0 )
    116232      {
    117233         MASS[i] = 52.08e0;
     234         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     235         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
     236         for( j = 0; j <= NLEV-1; j++ )
     237            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    118238      }
    119239      if( strcmp(corps[i], "C4H5") == 0 )
    120240      {
    121241         MASS[i] = 53.07e0;
     242         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     243         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
     244         for( j = 0; j <= NLEV-1; j++ )
     245            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    122246      }
    123247      if( strcmp(corps[i], "C4H6") == 0 )
    124248      {
    125249         MASS[i] = 54.09e0;
     250         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     251         sig     = 1.0e-16 * pow( ( siga + 4.6e0 ), 2.0e0 );
     252         for( j = 0; j <= NLEV-1; j++ )
     253            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    126254      }
    127255      if( strcmp(corps[i], "C4H10") == 0 )
    128256      {
    129257         MASS[i] = 58.13e0;
     258         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     259         sig     = 1.0e-16 * pow( ( siga + 4.687e0 ), 2.0e0 );
     260         for( j = 0; j <= NLEV-1; j++ )
     261            MD[i][j] = sqrt( p * TEMP[j] * m )
     262           / ( CT[j] * sig * omega(TEMP[j],epsa,531.4e0) );
    130263      }
    131264      if( strcmp(corps[i], "C6H") == 0 )
    132265      {
    133266         MASS[i] = 73.07e0;
     267         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     268         sig     = 1.0e-16 * pow( ( siga + 5.2e0 ), 2.0e0 );
     269         for( j = 0; j <= NLEV-1; j++ )
     270            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    134271      }
    135272      if( strcmp(corps[i], "C6H2") == 0 )
    136273      {
    137274         MASS[i] = 74.08e0;
     275         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     276         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
     277         for( j = 0; j <= NLEV-1; j++ )
     278            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    138279      }
    139280      if( strcmp(corps[i], "C8H2") == 0 )
    140281      {
    141282         MASS[i] = 98.10e0;
     283         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     284         sig     = 1.0e-16 * pow( ( siga + 6.0e0 ), 2.0e0 );
     285         for( j = 0; j <= NLEV-1; j++ )
     286            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    142287      }
    143288      if( strcmp( corps[i], "AC6H6" ) == 0 )
    144289      {
    145290         MASS[i] = 78.1136e0;
     291         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     292         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
     293         for( j = 0; j <= NLEV-1; j++ )   /* P. G. L. */
     294            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    146295      }
    147296      if( ( strcmp( corps[i], "C6H5" ) == 0 ) || ( strcmp( corps[i], "AC6H5" ) == 0 ) )
    148297      {
    149298         MASS[i] = 77.1136e0;
     299         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     300         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
     301         for( j = 0; j <= NLEV-1; j++ )
     302            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    150303      }
    151304      if( strcmp( corps[i], "C6H6" ) == 0 )
    152305      {
    153306         MASS[i] = 78.1136e0;
     307         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     308         sig     = 1.0e-16 * pow( ( siga + 5.4e0 ), 2.0e0 );
     309         for( j = 0; j <= NLEV-1; j++ )
     310            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    154311      }
    155312      if( strcmp(corps[i], "N2") == 0 )
    … …  
    160317      {
    161318         MASS[i] = 14.01e0;
     319         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     320         sig     = 1.0e-16 * pow( ( siga + 1.5e0 ), 2.0e0 );
     321         for( j = 0; j <= NLEV-1; j++ )
     322            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    162323      }
    163324      if( strcmp(corps[i], "NH") == 0 )
    164325      {
    165326         MASS[i] = 15.01e0;
     327         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     328         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
     329         for( j = 0; j <= NLEV-1; j++ )
     330            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    166331      }
    167332      if( strcmp(corps[i], "CN") == 0 )
    168333      {
    169334         MASS[i] = 26.02e0;
     335         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     336         sig     = 1.0e-16 * pow( ( siga + 3.2e0 ), 2.0e0 );
     337         for( j = 0; j <= NLEV-1; j++ )
     338            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    170339      }
    171340      if( strcmp(corps[i], "HCN") == 0 )
    172341      {
    173342         MASS[i] = 27.04e0;
     343         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     344         sig     = 1.0e-16 * pow( ( siga + 3.63e0 ), 2.0e0 );
     345         for( j = 0; j <= NLEV-1; j++ )
     346            MD[i][j] = sqrt( p * TEMP[j] * m )
     347           / ( CT[j] * sig * omega(TEMP[j],epsa,569.1e0) );
    174348      }
    175349      if( strcmp(corps[i], "H2CN") == 0 )
    176350      {
    177351         MASS[i] = 28.05e0;
     352         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     353         sig     = 1.0e-16 * pow( ( siga + 3.8e0 ), 2.0e0 );
     354         for( j = 0; j <= NLEV-1; j++ )
     355            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    178356      }
    179357      if( strcmp(corps[i], "C2N") == 0 )         /* C2N */
    180358      {
    181359         MASS[i] = 39.05e0;
     360         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     361         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     362         for( j = 0; j <= NLEV-1; j++ )
     363            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    182364      }
    183365      if( strcmp( corps[i], "CHCN" ) == 0 )
    184366      {
    185367         MASS[i]   = 39.05e0;
     368         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     369         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     370         for( j = 0; j <= NLEV-1; j++ )
     371            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    186372      }
    187373      if( strcmp( corps[i], "CH2CN" ) == 0 )
    188374      {
    189375         MASS[i]   = 40.04e0;
     376         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     377         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     378         for( j = 0; j <= NLEV-1; j++ )
     379            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    190380      }
    191381      if( strcmp( corps[i], "CH3CN" ) == 0 )
    192382      {
    193383         MASS[i]   = 41.05e0;
     384         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     385         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     386         for( j = 0; j <= NLEV-1; j++ )
     387            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    194388      }
    195389      if( strcmp( corps[i], "C2H3CN" ) == 0 )
    196390      {
    197391         MASS[i]   = 53.06e0;
     392         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     393         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     394         for( j = 0; j <= NLEV-1; j++ )
     395            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    198396      }
    199397      if( strcmp(corps[i], "NCCN") == 0 )        /* NCCN */
    200398      {
    201399         MASS[i] = 52.04e0;
     400         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     401         sig     = 1.0e-16 * pow( ( siga + 4.361e0 ), 2.0e0 );
     402         for( j = 0; j <= NLEV-1; j++ )
     403            MD[i][j] = sqrt( p * TEMP[j] * m )
     404           / ( CT[j] * sig * omega(TEMP[j],epsa,348.6e0) );
    202405      }
    203406      if( strcmp(corps[i], "C3N") == 0 )         /* C3N */
    204407      {
    205408         MASS[i] = 50.04e0;
     409         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     410         sig     = 1.0e-16 * pow( ( siga + 4.4e0 ), 2.0e0 );
     411         for( j = 0; j <= NLEV-1; j++ )
     412            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    206413      }
    207414      if( strcmp(corps[i], "HC3N") == 0 )        /* HC3N */
    208415      {
    209416         MASS[i] = 51.05e0;
     417         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     418         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
     419         for( j = 0; j <= NLEV-1; j++ )
     420            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    210421      }
    211422      if( strcmp( corps[i], "C4N2" ) == 0 )
    212423      {
    213424         MASS[i]   = 76.1e0;
     425         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     426         sig     = 1.0e-16 * pow( ( siga + 4.0e0 ), 2.0e0 );
     427         for( j = 0; j <= NLEV-1; j++ )
     428            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    214429      }
    215430      if( strcmp(corps[i], "H2O") == 0 )       
    216431      {
    217432         MASS[i] = 18.02e0;
     433         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     434         sig     = 1.0e-16 * pow( ( siga + 2.641e0 ), 2.0e0 );
     435         for( j = 0; j <= NLEV-1; j++ )
     436            MD[i][j] = sqrt( p * TEMP[j] * m )
     437           / ( CT[j] * sig * omega(TEMP[j],epsa,809.1e0) );
    218438      }
    219439      if( ( strcmp(corps[i], "O3P") == 0 ) || ( strcmp(corps[i], "O1D") == 0 ) )     
    220440      {
    221441         MASS[i] = 16.0e0;
     442         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     443         sig     = 1.0e-16 * pow( ( siga + 1.4e0 ), 2.0e0 );
     444         for( j = 0; j <= NLEV-1; j++ )
     445            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    222446      }
    223447      if( strcmp(corps[i], "OH") == 0 )       
    224448      {
    225449         MASS[i] = 17.01e0;
     450         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     451         sig     = 1.0e-16 * pow( ( siga + 3.0e0 ), 2.0e0 );
     452         for( j = 0; j <= NLEV-1; j++ )
     453            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
     454      }
     455      if( strcmp(corps[i], "HO2") == 0 )
     456      {
     457         MASS[i] = 33.01e0;
     458         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     459         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
     460         for( j = 0; j <= NLEV-1; j++ )
     461            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
     462      }
     463      if( strcmp(corps[i], "H2O2") == 0 )
     464      {
     465         MASS[i] = 33.01e0;
     466         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     467         sig     = 1.0e-16 * pow( ( siga + 3.5e0 ), 2.0e0 );
     468         for( j = 0; j <= NLEV-1; j++ )
     469            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
     470      }
     471      if( strcmp(corps[i], "O2") == 0 )
     472      {
     473         MASS[i] = 32.0e0;
     474         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     475         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
     476         for( j = 0; j <= NLEV-1; j++ )
     477            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
     478      }
     479      if( strcmp(corps[i], "O3") == 0 )
     480      {
     481         MASS[i] = 32.0e0;
     482         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     483         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
     484         for( j = 0; j <= NLEV-1; j++ )
     485            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    226486      }
    227487      if( strcmp(corps[i], "CO") == 0 )         
    228488      {
    229489         MASS[i] = 28.01e0;
     490         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     491         sig     = 1.0e-16 * pow( ( siga + 3.69e0 ), 2.0e0 );
     492         for( j = 0; j <= NLEV-1; j++ )
     493            MD[i][j] = sqrt( p * TEMP[j] * m )
     494           / ( CT[j] * sig * omega(TEMP[j],epsa,91.7e0) );
    230495      }
    231496      if( strcmp(corps[i], "HCO") == 0 )       
    232497      {
    233498         MASS[i] = 29.02e0;
     499         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     500         sig     = 1.0e-16 * pow( ( siga + 3.7e0 ), 2.0e0 );
     501         for( j = 0; j <= NLEV-1; j++ )
     502            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    234503      }
    235504      if( strcmp(corps[i], "CO2") == 0 )       
    236505      {
    237506         MASS[i] = 44.01e0;
     507         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     508         sig     = 1.0e-16 * pow( ( siga + 3.941e0 ), 2.0e0 );
     509         for( j = 0; j <= NLEV-1; j++ )
     510            MD[i][j] = sqrt( p * TEMP[j] * m )
     511           / ( CT[j] * sig * omega(TEMP[j],epsa,195.2e0) );
    238512      }
    239513      if( strcmp(corps[i], "CH2CO") == 0 )   
    240514      {
    241515         MASS[i] = 42.04e0;
     516         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     517         sig     = 1.0e-16 * pow( ( siga + 4.5e0 ), 2.0e0 );
     518         for( j = 0; j <= NLEV-1; j++ )
     519            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    242520      }
    243521      if( strcmp(corps[i], "CH2O") == 0 )   
    244522      {
    245523         MASS[i] = 30.03e0;
     524         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     525         sig     = 1.0e-16 * pow( ( siga + 3.75e0 ), 2.0e0 );
     526         for( j = 0; j <= NLEV-1; j++ )
     527            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    246528      }
    247529      if( ( strcmp(corps[i], "CH2OH") == 0 ) || ( strcmp(corps[i], "CH3O") == 0 ) ) 
    248530      {
    249531         MASS[i] = 31.04e0;
     532         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     533         sig     = 1.0e-16 * pow( ( siga + 3.4e0 ), 2.0e0 );
     534         for( j = 0; j <= NLEV-1; j++ )
     535            MD[i][j] = sqrt( p * TEMP[j] * m ) / ( CT[j] * sig );
    250536      }
    251537      if( strcmp(corps[i], "CH3OH") == 0 )       
    252538      {
    253539         MASS[i] = 32.042e0;
     540         m       = ( ma + MASS[i] ) / ( ma * MASS[i] );
     541         sig     = 1.0e-16 * pow( ( siga + 3.626e0 ), 2.0e0 );
     542         for( j = 0; j <= NLEV-1; j++ )
     543            MD[i][j] = sqrt( p * TEMP[j] * m )
     544           / ( CT[j] * sig * omega(TEMP[j],epsa,481.8e0) );
    254545      }
    255546   }

  • trunk/LMDZ.TITAN/libf/chimtitan/disso.c

    r3 r1126  
    77#include "titan.h"
    88
    9 void disso_( double KRPD[][NLEV][RDISS+1][15], int *NLAT )
     9void disso_( double KRPD[][NLRT][RDISS+1][15], int *NLAT )
    1010{
    1111   static double sH2[62] = {  /* incertain en dessous de 70 et en dessus de 85... */
    … …  
    179179   double f,**flux;
    180180   double flact;
    181    char   name[60];
     181   char   name[60],dir[24];
    182182   FILE   *fp,*out;
    183183
    184    flux = dm2d(0,NLEV,0,14);
    185 
     184   flux = dm2d(0,NLRT,0,14);
     185
     186/* lecture des flux actiniques:
     187   - suppose que l'executable est dans $LMDGCM/RUN/xxx/
     188   - et les moyennes dans $LMDGCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
     189*/
     190   strcpy( dir, "../../INPUT/PHOT" );
     191   if( (*NLAT) < 10 )
     192   {
     193     strcat( dir, "0"  );
     194     strcat( dir, (const char *)ecvt((float)(*NLAT),1,&x,&x) );
     195   }
     196   else
     197     strcat( dir, (const char *)ecvt((float)(*NLAT),2,&x,&x) );
     198   strcat( dir, "x/Moy_"  );
     199   printf( "Directories for actinic fluxes: %s \n", dir );
     200 
    186201   for( lat = 0; lat <= (*NLAT)-1; lat++ )   /* Old array is set equal to 0. */
    187      for( j = 0; j <= NLEV-1; j++ )
     202     for( j = 0; j <= NLRT-1; j++ )
    188203       for( i = 0; i <= RDISS; i++ )                 
    189204         for( s = 0; s <= 14; s++ )
    … …  
    195210   for( i = 0; i <= 16; i++ ) sC2H6[i] = 0.0e0;
    196211
    197 //   for( lat = 25; lat <= 25; lat++ )     /* Main loop on latitude */
    198212   for( lat = 0; lat <= (*NLAT)-1; lat++ ) /*     Main loop on latitude */
    199213   for( i = 10; i <= 310; i += 5 )             /* Main loop on wavelength. */
    200214   {
    201 /* lecture des flux actiniques:
    202    - suppose que l'executable est dans $LMDGCM/RUN/xxx/
    203    - et les moyennes dans GCCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
    204    (sauf a l'idris...)
    205 */
    206       strcpy( name, "../../INPUT/PHOT49/Moy_" );
    207 //    strcpy( name, "PHOT49/Moy_" );
     215      strcpy( name, dir );
    208216      if( (lat+1) < 10 )
    209217      {
    … …  
    213221      else
    214222       strcat( name, (const char *)ecvt((float)(lat+1),2,&x,&x) );
    215       if( i < 160 ) strcat( name, "/photgcm3a" );
    216       else          strcat( name, "/photgcm3a" );
     223      if( i < 160 ) strcat( name, "/photmoy3a" );
     224      else          strcat( name, "/photmoy3a" );
    217225      if( i < 100 )
    218226      {
    … …  
    232240         exit(0);
    233241      }
    234       for( j = 0; j <= NLEV-1; j++ )
     242      for( j = 0; j <= NLRT-1; j++ )
    235243      {
    236244         fscanf( fp,"%d ",&l );
    … …  
    248256
    249257      for( s = 0; s <= 14; s++ )
    250        for( j = 0; j <= NLEV-1; j++ )
     258       for( j = 0; j <= NLRT-1; j++ )
    251259       {
    252260         f = flux[j][s] * sol[l] / ( 9.5e0 * 9.5e0 );   /* !! # de reac de 0 a RDISS-1 !! */
    … …  
    369377    for( s = 0; s <= 14; s++ )
    370378    {
    371      for( j = 36; j <= NLEV-1; j++ )  /* level 36 ~ 200 km */
     379     for( j = 99; j <= NLRT-1; j++ )  /* level 100 = 200 km */
    372380      KRPD[lat][j][RDISS][s] = 1.0e-16;
    373      for( j = 26; j <= 35; j++ )      /* level 26 ~ 100 km */
    374       KRPD[lat][j][RDISS][s] = 1.0e-17+9.e-18*(j-26);
    375      for( j = 23; j <= 25; j++ )      /* level 23 ~ 70 km */
    376       KRPD[lat][j][RDISS][s] = pow(10.,(-23+2*(j-23)));
    377      for( j = 0; j <= 22; j++ )
     381     for( j = 49; j <= 98; j++ )      /* level 50 = 100 km */
     382      KRPD[lat][j][RDISS][s] = 1.0e-17+1.8e-18*(j-49);
     383     for( j = 34; j <= 48; j++ )      /* level 35 = 70 km */
     384      KRPD[lat][j][RDISS][s] = pow(10.,(-23+0.4*(j-34)));
     385     for( j = 0; j <= 33; j++ )
    378386      KRPD[lat][j][RDISS][s] = 0.0e0;
    379387    }

  • trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c

    r1057 r1126  
    1313void gptitan_(
    1414     double *RA, double *TEMP, double *NB,
    15      char CORPS[][10], double Y[][NLEV], double *FTOP,
    16      double *DECLIN, double *FIN, int *LAT, double *MASS,
    17      double *botCH4,
    18      double KRPD[][NLEV][RDISS+1][15], double KRATE[][NLEV],
     15     char CORPS[][10], double Y[][NLEV],
     16     double *FIN, int *LAT, double *MASS, double MD[][NLEV],
     17     double *KEDD, double *botCH4, double KRATE[][NLEV],
    1918     int reactif[][5], int *nom_prod, int *nom_perte,
    2019     int prod[][200], int perte[][200][2], int *aerprod, int *utilaer,
    … …  
    2423{
    2524   char   outlog[100],corps[100][10];
    26    int    dec,declinint,i,j,k,l;
     25   int    i,j,k,l;
    2726   int    ireac,ncom1,ncom2;
    28    double  *fl,*fp,*mu,**jac,*ym1;
    29    double  *fluxtop,fluxCH4;
    30    double  cm,conv,cp,delta,deltamax,dv,dr,drp,drm;
    31    double  rr,np,nm,factdec,s,test,time,ts,v;
    32    double *fd,**jacd;
     27   double  ***a,***b,**c;
     28   double  *fl,*fp,*mu,**jac,**ym1,**f;
     29   double  fluxCH4;
     30   double  conv,delta,deltamax;
     31   double  cm,cp,dim,dip,dm,dp,dym,dyp,km,kp,r,dra,dram,drap;
     32   double  np,nm,s,test,time,ts,v,dv;
    3333   char   str2[15];
    3434   FILE   *out;
    … …  
    5353
    5454/* DEBUG */
    55       printf("CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
     55      printf("CHIMIE: lat=%d\n",(*LAT)+1);
    5656/**/
    5757
    … …  
    6464   strcpy( outlog, "chimietitan" );
    6565   strcat( outlog, ".log" );
     66   out = fopen( outlog, "w" );
     67   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
     68   fclose( out );
     69
    6670   deltamax = 1.e5;
     71   test = 1.0e-15;
     72
     73/* valeur de r:
     74r = g0 R0^2 / R * 2 * 1E-3
     75avec g0 en cm/s2, R0 en km, mu et mass en g
     76*/
     77   r        = 21.595656e0;
    6778
    6879/*  DEBUG
    6980            out = fopen( outlog, "a" );
    70    fprintf(out,"CHIMIE: lat=%d declin=%e\n",(*LAT)+1,(*DECLIN));
     81   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
    7182            fclose( out );
    7283*/
    73    ym1     = dm1d( 0,   NC-1 );
    7484   fl      = dm1d( 0,   NC-1 );
    7585   fp      = dm1d( 0,   NC-1 );
    76    fd      = dm1d( 0,   NC-1 );
    7786   mu      = dm1d( 0, NLEV-1 );
    78    fluxtop = dm1d( 0,   NC-1 );
     87   ym1     = dm2d( 0,   NC-1, 0, NLEV-1 );
     88   f       = dm2d( 0,   NC-1, 0, NLEV-1 );
    7989   jac     = dm2d( 0,   NC-1, 0,   NC-1 );
    80    jacd    = dm2d( 0,   NC-1, 0,   NC-1 );
     90   c       = dm2d( 0, NLEV-1, 0,   NC-1 );
     91   a       = dm3d( 0, NLEV-1, 0,   NC-1, 0,   NC-1 );
     92   b       = dm3d( 0, NLEV-1, 0,   NC-1, 1,    2 );
    8193
    8294/* DEBUG */
    … …  
    8799*/
    88100
    89 /* initialisation krate pour dissociations */
    90 
    91    if( ( (*DECLIN) *10 + 267 ) < 14. )
    92    {
    93      declinint = 0;
    94      dec = 0;
    95    }
    96    else
    97    { 
    98        if( ( (*DECLIN) * 10 + 267 ) > 520. )
    99        {
    100           declinint = 14;
    101           dec = 534;
    102        }
    103        else
    104        {
    105           declinint = 1;
    106           dec       = 27;
    107           while( ( (*DECLIN) * 10 + 267 ) >= (float)(dec+20) )
    108           {
    109             dec += 40;
    110             declinint++;
    111           }
    112        }
    113    }
    114    if( ( (*DECLIN) >= -24. ) && ( (*DECLIN) <= 24. ) )
    115            factdec = ( (*DECLIN) - (dec-267)/10. ) / 4.;
    116    else
    117            factdec = ( (*DECLIN) - (dec-267)/10. ) / 2.7;
    118 
    119    for( i = 0; i <= RDISS; i++ )   /* RDISS correspond a la dissociation de N2 par les GCR */
    120       for( j = 0; j <= NLEV-1; j++ )
    121       {
    122          if( factdec < 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint-1] * fabs(factdec)
    123                                          + KRPD[*LAT][j][i][declinint] * ( 1 + factdec);
    124          if( factdec > 0. )  KRATE[i][j] = KRPD[*LAT][j][i][declinint+1] * factdec
    125                                          + KRPD[*LAT][j][i][declinint] * ( 1 - factdec);
    126          if( factdec == 0. ) KRATE[i][j] = KRPD[*LAT][j][i][declinint];
    127       }
    128        
    129101/* initialisation mu, CH4 au sol */
    130102     
    … …  
    143115   }
    144116
     117/* initialisation compo avant calcul */
     118   for( j = NLEV-1; j >= 0; j-- )
     119      for( i = 0; i <= ST-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
     120
     121/*
     122==========================================================================
     123   STRATEGIE:
     124    INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
     125    PUIS INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
     126==========================================================================
     127
     128PREMIERE ETAPE:
     129===============
     130    INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
     131===============
     132*/
     133
    145134/* ****************** */
    146 /*  Main loop: level  */
     135/*  Time loop:        */
    147136/* ****************** */
    148137
    149    for( j = NLEV-1; j >= 0; j-- )
     138time     = ts = 0.0e0;
     139delta    = 1.e-3;
     140
     141while( time < (*FIN) )     
     142{
     143
     144
     145/* DEBUG
     146   for( j = NLEV-1; j >= NLD; j-- )
    150147   {
    151 
    152 /* DEBUG
    153148            out = fopen( outlog, "a" );
    154149        fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
    … …  
    160155    for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
    161156            fclose( out );
    162 
    163     printf("%d %e %e %e\n",declinint,(RA[j]-R0),NB[j],TEMP[j]);
    164     for (i=0;i<=RDISS-1;i++) printf("%d %e\n",i,KRPD[*LAT][j][i][declinint]);
    165     for (i=0;i<=ST-1;i++)    printf("%10s %e\n",corps[i],FTOP[i]);
    166 
     157    }
    167158   exit(0);
    168159*/
    169160
    170       time     = ts = 0.0e0;
    171 /*      delta    = (*FIN);  */
    172       delta    = 1.e-3;
    173 
    174       for( i = 0; i <= ST-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
    175 
    176 /* ++++++++++++ */
    177 /*  time loop.  */
    178 /* ++++++++++++ */
    179 
    180       while( time < (*FIN) )     
    181       {
    182 
    183 /* Calcul de f et de la matrice jacobienne */
    184 /* --------------------------------------- */
    185 
    186          for( i = 0; i <= ST-1; i++ )     /* productions et pertes chimiques */
     161
     162/* ------------------------------ */
     163/* Calculs variations et jacobien */
     164/* ------------------------------ */
     165
     166   for( j = NLEV-1; j >= NLD; j-- )
     167   {
     168
     169/* init of step */
     170/* ------------ */
     171         for( i = 0; i <= ST-1; i++ )
     172         {
     173            fp[i] = fl[i] = 0.0e0;
     174            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
     175         }
     176
     177/* Chimie */
     178/* ------ */
     179
     180/* productions et pertes chimiques */
     181         for( i = 0; i <= ST-1; i++ )   
    187182         {
    188183            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
    189 
    190             fp[i] = fl[i] = 0.0e0;                    /* init for next step */
    191             for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
    192184
    193185            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
    … …  
    226218            }
    227219         }
     220
    228221
    229222/* Aerosols */
    … …  
    283276     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
    284277         }
    285            
     278     
     279       
    286280/* H -> H2 on haze particles */
    287281/* ------------------------- */
    … …  
    299293*/
    300294
    301               fp[utilaer[0]] += ( dyh  * NB[j] );
     295              fp[utilaer[0]] += ( dyh  * NB[j] );
     296   /* pourquoi pas *2 ?? cf gptit dans 2da... */
     297
    302298              fp[utilaer[1]] += ( dyh2 * NB[j] );
    303299              if( Y[utilaer[0]][j] != 0.0 )
    304300       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
    305          }
    306 
    307          for( i = 0; i <= ST-1; i++ )     /* finition pour inversion */
    308          {
     301   /* pourquoi pas *2 ?? cf gptit dans 2da... */
     302         }
     303
     304
     305/* Backup jacobian level j. */
     306/* ------------------------ */
     307         for( i = 0; i <= ST-1; i++ )
    309308            for( k = 0; k <= ST-1; k++ )
    310             {
    311                jac[i][k] *= ( -THETA * delta );     /* Correction time step. */
    312                if( k == i ) jac[k][k] += NB[j];     /* Correction diagonal. */
    313                jacd[i][k] = (double)jac[i][k];
    314             }
    315            
    316             fd[i] = (double)(delta * ( fp[i] - Y[i][j]*fl[i] ));
    317          }
    318 
    319 /*         for( i = ST; i <= NC-1; i++ )    pas d'inversion (soot,prod): que faire?
    320          {
    321             Y[i][j] = ??? ;
    322          }
    323 */
    324 
     309               a[j][i][k] = jac[i][k];     
     310
     311
     312/* Diffusion verticale et flux exterieurs */
     313/* -------------------------------------- */
     314
     315/*
     316pour dy/dr, dr doit etre en cm...
     317pareil pour dphi/dr
     318*/
     319         for( i = 0; i <= ST-1; i++ )
     320         {
     321
     322/* First level. */
     323            if( j == NLD )
     324            {
     325               v = dv = 0.0e0;
     326               dra = RA[j+1]-RA[j];
     327
     328               cp  = (NB[j+1]+NB[j])/2.;  /* Mean total concentration. */
     329               dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
     330                     pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
     331               dp  = (MD[i][j]+MD[i][j+1])/2.;     /* Mean molecular diffusion. */
     332               dyp = (Y[i][j+1]-Y[i][j])/(RA[j+2]-RA[j])*2.e-5; /* Delta y level +1. */
     333               kp  = (KEDD[j+1]+KEDD[j])/2.;       /* Mean eddy diffusion. */
     334             /* div phi. */
     335               f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
     336                              + kp * dyp )
     337                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.))
     338                       + fp[i] - Y[i][j]*fl[i] + v;
     339             /* dphi / dy this level. */
     340               a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
     341                                    - 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
     342                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) + dv );
     343             /* dphi / dy level +1. */
     344               c[j][i]     = -THETA * delta
     345                             * cp * ( dp * 0.5e0 * dip
     346                                    + 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
     347                        * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.));
     348            }
     349/* Last level. */
     350            else if( j == NLEV-1 )
     351            {
     352               v = dv = 0.0e0;
     353               dra = RA[NLEV-1]-RA[NLEV-2];
     354
     355   /* Jeans escape */
     356               if( strcmp(corps[i], "H") == 0 ) 
     357               {
     358                 dv = top_H  * NB[NLEV-1]
     359                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
     360                 v  = dv * Y[i][NLEV-1];
     361               }
     362               if( strcmp(corps[i], "H2") == 0 )
     363               {
     364                 dv = top_H2 * NB[NLEV-1]
     365                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
     366                 v  = dv * Y[i][NLEV-1];
     367               }
     368   /* Input flux for N(4S) */
     369               if( strcmp(corps[i], "N4S") == 0 )
     370                 v  = top_N4S
     371                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
     372
     373               cm  = (NB[NLEV-1]+NB[NLEV-2])/2.;  /* Mean total concentration. */
     374               dim = r * (MASS[i]-(mu[NLEV-1]+mu[NLEV-2])/2.)
     375                       / (TEMP[NLEV-1]+TEMP[NLEV-2])
     376                       / pow( RA[NLEV-1],   2.0e0 );  /* Delta i,j level -1. */
     377               dm  = (MD[i][NLEV-1]+MD[i][NLEV-2])/2.;    /* Mean molecular diffusion. */
     378               dym = (Y[i][NLEV-1]-Y[i][NLEV-2])/dra*1.e-5; /* Delta y level -1. */
     379               km  = (KEDD[NLEV-1]+KEDD[NLEV-2])/2.;      /* Mean eddy diffusion. */
     380             /* div phi. */
     381               f[i][NLEV-1] = fp[i] - Y[i][NLEV-1]*fl[i] - v
     382                       - cm * ( dm * ( (Y[i][NLEV-1]+Y[i][NLEV-2])/2. * dim + dym )
     383                              + km * dym )
     384                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
     385             /* dphi / dy this level */
     386               a[NLEV-1][i][i] -= ( cm * ( dm * 0.5e0 * dim
     387                                    + 1.e-5/dra * (dm + km ) )
     388                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)) + dv );
     389             /* dphi / dy level -1. */
     390               b[NLEV-1][i][2]  =  THETA * delta
     391                                  * cm * ( dm * 0.5e0 * dim
     392                                    - 1.e-5/dra * (dm + km ) )
     393                        * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
     394            }
     395            else
     396            {
     397               v = dv = 0.0e0;
     398               dram=(RA[j+1]-RA[j-1])/2.;
     399               if (j<NLEV-2)
     400                 drap=(RA[j+1]-RA[j-1])/2.;
     401               else
     402                 drap=dram;
     403
     404               cm  = (NB[j]+NB[j-1])/2.;       /* Mean concentration level -1. */
     405               cp  = (NB[j]+NB[j+1])/2.;       /* Mean concentration level +1. */
     406               dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
     407                     pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
     408               dim = r * (MASS[i]-(mu[j]+mu[j-1])/2.) / (TEMP[j]+TEMP[j-1]) /
     409                     pow( RA[j],   2.0e0 );    /* Delta i,j level -1. */
     410               dm  = (MD[i][j-1]+MD[i][j])/2.;    /* Mean molecular diffusion level -1. */
     411               dp  = (MD[i][j+1]+MD[i][j])/2.;    /* Mean molecular diffusion level +1. */
     412               dym = (Y[i][j]-Y[i][j-1])/dram*1.e-5; /* Delta y level -1. */
     413               dyp = (Y[i][j+1]-Y[i][j])/drap*1.e-5; /* Delta y level +1. */
     414               km  = (KEDD[j]+KEDD[j-1])/2.;      /* Mean eddy diffusion level -1. */
     415               kp  = (KEDD[j]+KEDD[j+1])/2.;      /* Mean eddy diffusion level +1. */
     416             /* div phi. */
     417               f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
     418                              + kp * dyp )
     419                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
     420                       - cm * ( dm * ( (Y[i][j]+Y[i][j-1])/2. * dim + dym )
     421                              + km * dym )
     422                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.))
     423                       + fp[i] - fl[i] * Y[i][j] + v;
     424             /* dphi / dy this level */
     425               a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
     426                                    - 1.e-5/drap * (dp + kp) )
     427                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
     428                             - cm * ( dm * 0.5e0 * dim
     429                                    + 1.e-5/dram * (dm + km ) )
     430                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) );
     431             /* dphi / dy level -1. */
     432               b[j][i][2]  =  THETA * delta
     433                             * cm * ( dm * 0.5e0 * dim
     434                                    - 1.e-5/dram * (dm + km ) )
     435                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.));
     436             /* dphi / dy level +1. */
     437               c[j][i]     = -THETA * delta
     438                             * cp * ( dp * 0.5e0 * dip
     439                                    + 1.e-5/drap * (dp + kp) )
     440                        * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.));
     441            }
     442         }
     443
     444
     445
     446/* finition pour inversion */
     447/* ----------------------- */
     448
     449         for( i = 0; i <= ST-1; i++ )
     450         {
     451            for( k = 0; k <= ST-1; k++ )
     452            {
     453               a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
     454               if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
     455            }
     456            f[i][j] *= delta;
     457         }
     458
     459   }
     460
     461
     462/* -------------------------------- */
    325463/* Inversion of matrix cf method LU */
    326464/* -------------------------------- */
    327465
     466   for( j = NLD+1; j <= NLEV-1; j++ )
     467   {
     468         solve( a, j-1, 0, ST-1 );
     469         for( i = 0; i <= ST-1; i++ )
     470         {
     471            s = 0.0e0;
     472            for( k = 0; k <= ST-1; k++ )
     473            {
     474               a[j][i][k] -= ( b[j][i][2] * c[j-1][k] * a[j-1][i][k] );
     475               s          += ( b[j][i][2] * f[k][j-1] * a[j-1][i][k] );
     476            }
     477            f[i][j] -= s;
     478         }
     479   }
     480   solve( a, NLEV-1, 0, ST-1 );
     481   for( j = NLEV-1; j >= NLD; j-- )     
     482   {
     483         if( j != NLEV-1 )
     484            for( i = 0; i <= ST-1; i++ ) f[i][j] -= ( c[j][i] * b[j+1][i][1] );
     485         for( i = 0; i <= ST-1; i++ )
     486         {
     487            s = 0.0e0;
     488            for( k = 0; k <= ST-1; k++ ) s += ( a[j][i][k] * f[k][j] );
     489            b[j][i][1]  = s;
     490            Y[i][j]    += s;
     491            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
     492         }
     493   }
     494
     495/* ------------------ */
     496/* Tests et evolution */
     497/* ------------------ */
     498
     499/* Calcul deviation */
     500/* ---------------- */
     501
     502   for( j = NLD; j <= NLEV-1; j++ )
     503      for( i = 0; i <= ST-1; i++ )
     504         if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
     505         {
     506               conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
     507               if( conv > ts )
     508               {
     509/*
     510                  if( conv >= 0.1 )
     511                  {
     512                     out = fopen( outlog, "a" );
     513                     fprintf( out, "Lat no %d;", (*LAT)+1);
     514                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
     515                     fclose( out );
     516                  }
     517*/
     518                  ts = conv;
     519               }
     520         }
     521
     522/* test deviation */
     523/* -------------- */
     524
     525         if( ts < 0.1e0 )
     526         {
     527            for( i = 0; i <= ST-1; i++ )
     528               for( j = NLD; j <= NLEV-1; j++ )
     529                 if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
     530                 {
     531                  out = fopen( outlog, "a" );
     532                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
     533                           corps[i], ym1[i], Y[i][j], j );
     534                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
     535                  fclose( out );
     536                  exit(0);
     537//                  Y[i][j] = 1.e-20;
     538                 }
     539            for( j = NLD; j <= NLEV-1; j++ )
     540               for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
     541            time += delta;
     542            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
     543            if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
     544            if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
     545            if( ( ts > 1.00e-3 ) && ( ts < 5.0e-3 ) ) delta *= 3.0e0;
     546            if( ( ts > 5.00e-3 ) && ( ts < 0.01e0 ) ) delta *= 1.5e0;
     547            if( ( ts > 0.010e0 ) && ( ts < 0.03e0 ) ) delta *= 1.2e0;
     548            if( ( ts > 0.030e0 ) && ( ts < 0.05e0 ) ) delta *= 1.1e0;
     549
     550//            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
     551//            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
     552         
     553            delta = min( deltamax, delta );
     554         }
     555         else
     556         {
     557            for( j = NLD; j <= NLEV-1; j++ )
     558               for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
     559
     560            if(   ts > 0.8 )                    delta *= 1.e-6;
     561            if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
     562            if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
     563            if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
     564            if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
     565            if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
     566         }
     567         ts = 0.0e0;
     568
     569         out = fopen( outlog, "a" );
     570         fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
     571         fclose( out );
     572
     573}
     574/* **************** */       
     575/* end of time loop */
     576/* **************** */       
     577
     578/*
     579==========================================================================
     580
     581SECONDE ETAPE:
     582===============
     583    INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
     584===============
     585*/
     586   if( NLD != 0 )
     587   for( j = NLD-1; j >= 0; j-- )
     588   {
     589      time     = ts = 0.0e0;
     590      delta    = 1.e-3;
     591
     592/* ++++++++++++ */
     593/*  time loop.  */
     594/* ++++++++++++ */
     595
     596      while( time < (*FIN) )     
     597      {
     598
     599/* init of step */
     600/* ------------ */
     601         for( i = 0; i <= ST-1; i++ )
     602         {
     603            fp[i] = fl[i] = 0.0e0;
     604            for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
     605         }
     606
     607/* Chimie */
     608/* ------ */
     609
     610/* productions et pertes chimiques */
     611         for( i = 0; i <= ST-1; i++ )   
     612         {
     613            Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
     614
     615            for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
     616            {
     617               ireac = prod[i][l];                  /* Number of the reaction involves. */
     618               ncom1 = reactif[ireac][0];           /* First compound which reacts. */
     619               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     620               {
     621                  jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
     622                  fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
     623               }
     624               else                                 /* General case. */
     625               {
     626                  ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
     627                  jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
     628                  jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
     629                  fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
     630               }
     631            }
     632           
     633            for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
     634            {
     635               ireac = perte[i][l][0];              /* Reaction number. */
     636               ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
     637               if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     638               {
     639                  jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
     640                  fl[i]     += ( KRATE[ireac][j] * NB[j] );
     641               }
     642               else                                 /* General case. */
     643               {
     644                  jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
     645                  jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
     646                  fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
     647               }
     648            }
     649         }
     650
     651
     652/* Aerosols */
     653/* -------- */
     654         if( (*aerprod) == 1 )
     655         {
     656             aer(corps,TEMP,NB,Y,&j,k_dep,faer,
     657              &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
     658              mmolaer,productaer,csurn,csurh);
     659
     660             for( i = 0; i <= 3; i++ )
     661             {
     662               PRODAER[i][j] = productaer[i];
     663                  MAER[i][j] = mmolaer[i];
     664                   CSN[i][j] = csurn[i];
     665                   CSH[i][j] = csurh[i];
     666             }
     667/* DEBUG
     668printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
     669if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
     670      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
     671              fp[utilaer[2]],dyc2h2*NB[j]);
     672if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
     673      printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
     674              fp[utilaer[5]],dyhcn*NB[j]);
     675if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
     676      printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
     677              fp[utilaer[6]],dyhc3n*NB[j]);
     678if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
     679      printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
     680              fp[utilaer[13]],dynccn*NB[j]);
     681if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
     682      printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
     683              fp[utilaer[14]],dych3cn*NB[j]);
     684if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
     685      printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
     686              fp[utilaer[15]],dyc2h3cn*NB[j]);
     687*/
     688
     689             fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
     690             fp[utilaer[5]] -= (    dyhcn * NB[j] );
     691             fp[utilaer[6]] -= (   dyhc3n * NB[j] );
     692             fp[utilaer[13]]-= (   dynccn * NB[j] );
     693             fp[utilaer[14]]-= (  dych3cn * NB[j] );
     694             fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
     695             if( Y[utilaer[2]][j]  != 0.0 )
     696       jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
     697             if( Y[utilaer[5]][j]  != 0.0 )
     698       jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
     699             if( Y[utilaer[6]][j]  != 0.0 )
     700       jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
     701             if( Y[utilaer[13]][j] != 0.0 )
     702     jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
     703             if( Y[utilaer[14]][j] != 0.0 )
     704     jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
     705             if( Y[utilaer[15]][j] != 0.0 )
     706     jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
     707         }
     708     
     709       
     710/* H -> H2 on haze particles */
     711/* ------------------------- */
     712         if( (*htoh2) == 1 )
     713         {
     714              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
     715/* dyh <= 0 / 1.0 en adsor., 1 en reac. */
     716
     717/* DEBUG
     718printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
     719if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
     720printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
     721if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
     722printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
     723*/
     724
     725              fp[utilaer[0]] += ( dyh  * NB[j] );
     726   /* pourquoi pas *2 ?? cf gptit dans 2da... */
     727
     728              fp[utilaer[1]] += ( dyh2 * NB[j] );
     729              if( Y[utilaer[0]][j] != 0.0 )
     730       jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
     731   /* pourquoi pas *2 ?? cf gptit dans 2da... */
     732         }
     733
     734
     735/* Backup jacobian level j. */
     736/* ------------------------ */
     737         for( i = 0; i <= ST-1; i++ )
     738         {
     739            for( k = 0; k <= ST-1; k++ )
     740               a[j][i][k] = jac[i][k];     
     741            f[i][j] = fp[i] - fl[i] * Y[i][j];
     742         }
     743
     744
     745/* finition pour inversion */
     746/* ----------------------- */
     747
     748         for( i = 0; i <= ST-1; i++ )
     749         {
     750            for( k = 0; k <= ST-1; k++ )
     751            {
     752               a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
     753               if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
     754            }
     755            f[i][j] *= delta;
     756         }
     757
     758
     759/* Inversion of matrix cf method LU */
     760/* -------------------------------- */
     761
    328762/* inversion by blocs: */
    329763/* Hydrocarbons */
    330764
    331          solve( jacd, fd, 0, NHC-1 );             
    332 
     765         solve_b( a, f, j, 0, NHC-1 );             
    333766         for( i = 0; i <= NHC-1; i++ )
    334767         {
    335             Y[i][j] += (float)fd[i];
     768            Y[i][j] += f[i][j];
    336769            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
    337770         }
    … …  
    339772/* Nitriles */
    340773
    341          solve( jacd, fd, NHC, ST-1 );             
    342 
     774         solve_b( a, f, j, NHC, ST-1 );             
    343775         for( i = NHC+1; i <= ST-1; i++ )
    344776         {
    345             Y[i][j] += (float)fd[i];
     777            Y[i][j] += f[i][j];
    346778            if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
    347779         }
    348780
    349781/* end inversion by blocs: */
    350 
    351          for( i = 0; i <= ST-1; i++ )
    352          {
    353782
    354783/* CH4 au sol */
    355784/* ---------- */
    356 
    357             if( ( strcmp(corps[i], "CH4") == 0 ) && ( j == 0 ) && ( Y[i][j] < *botCH4 ) )
    358             {
    359                 fluxCH4 += (*botCH4 - Y[i][0]);
    360                 Y[i][0] = *botCH4;
    361             }
    362 
    363          }
    364          
    365 /* test evolution delta */
    366 /* -------------------- */
     785   for( i = 0; i <= ST-1; i++ )
     786     if( ( strcmp(corps[i], "CH4") == 0 ) && (j==0) && ( Y[i][0] < *botCH4 ) )
     787     {
     788          fluxCH4 += (*botCH4 - Y[i][0]);
     789           Y[i][0] = *botCH4;
     790     }
     791
     792/* ------------------ */
     793/* Tests et evolution */
     794/* ------------------ */
     795
     796/* Calcul deviation */
     797/* ---------------- */
    367798
    368799         for( i = 0; i <= ST-1; i++ )
    369800         {
    370801            test = 1.0e-15;
    371             if( ( Y[i][j] > test ) && ( ym1[i] > test ) )
    372             {
    373                conv = fabs( Y[i][j] - ym1[i] ) / ym1[i];
     802            if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
     803            {
     804               conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
    374805
    375806               if( conv > ts )
    … …  
    389820         }
    390821
     822/* test deviation */
     823/* -------------- */
     824
    391825         if( ts < 0.1e0 )
    392826         {
    … …  
    396830                  out = fopen( outlog, "a" );
    397831                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
    398                            corps[i], ym1[i], Y[i][j], j );
    399                   for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i],Y[i][k] );
     832                           corps[i], ym1[i][j], Y[i][j], j );
     833                  for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i][j],Y[i][k] );
    400834                  fclose( out );
    401835             //     exit(0);
    402836                  Y[i][j] = 1.e-20;
    403837                 }
    404             for( i = 0; i <= NC-1; i++ ) ym1[i] = max(Y[i][j],1.e-30);
     838            for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
    405839            time += delta;
    406840            if(   ts < 1.00e-5 )                      delta *= 1.0e2;
    … …  
    414848         else
    415849         {
    416             for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i];
     850            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
    417851
    418852            if(   ts > 0.8 )                    delta *= 1.e-6;
    … …  
    439873            fluxCH4 *= ( MASS[i]/(6.022e23*time) );
    440874
    441    }
    442    
    443 /* **************** */       
    444 /* end of main loop */
    445 /* **************** */       
    446 
    447 /* Plafond: !! OU !!  flux vertical     */
    448 /* ------------------------------------ */
    449 
    450    for( i = 0; i <= ST-1; i++ )
    451       if( FTOP[i] != 0.0e0 )
    452       {
    453              fluxtop[i]    = (- FTOP[i]/NB[NLEV-2]) * MASS[i]/6.022e23;
    454              Y[i][NLEV-2] += FTOP[i]/NB[NLEV-2]*(*FIN);
    455              Y[i][NLEV-2]  = max(Y[i][NLEV-2],0.0e0);
    456 // on ajuste aussi le niveau dans la derniere couche...
    457 // pour eviter les effets vers le haut
    458              Y[i][NLEV-1]  = Y[i][NLEV-2];
    459       }
     875   }  /*  boucle j */
     876
     877
     878/*
     879==========================================================================
     880
     881FINALISATION:
     882===============
     883*/
     884      for( i = 0; i <= ST-1; i++ )
     885         if( strcmp(corps[i],"CH4") == 0 )
     886            fluxCH4 *= ( MASS[i]/(6.022e23*time) );
    460887
    461888/* Niveau de N2 */
    … …  
    465892   {
    466893      conv = 0.0e0;
    467       for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
    468       for( i = 0; i <= ST-1; i++ ) if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
     894      for( i = 0; i <= ST-1; i++ )
     895         if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
     896      for( i = 0; i <= ST-1; i++ )
     897         if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
    469898   }
    470899
    … …  
    479908   }
    480909
    481    fdm1d(     ym1, 0 );
    482910   fdm1d(      fl, 0 );
    483911   fdm1d(      fp, 0 );
    484    fdm1d(      fd, 0 );
    485912   fdm1d(      mu, 0 );
    486    fdm1d( fluxtop, 0 );
     913   fdm2d(     ym1, 0,   NC-1, 0 );
     914   fdm2d(       f, 0,   NC-1, 0 );
    487915   fdm2d(     jac, 0,   NC-1, 0 );
    488    fdm2d(    jacd, 0,   NC-1, 0 );
     916   fdm2d(       c, 0, NLEV-1, 0 );
     917   fdm3d(       a, 0, NLEV-1, 0,   NC-1, 0 );
     918   fdm3d(       b, 0, NLEV-1, 0,   NC-1, 1 );
    489919}

  • trunk/LMDZ.TITAN/libf/chimtitan/solve.c

    r3 r1126  
    22/* cf Numerical Recipes LU Method for equations numbers */
    33/* GCCM */
    4 /* une seule cellule concernee */
     4/* similaire a inv (GP) */
    55/* la matrice a est inversee seulement sur le bloc [n0;n1][n0;n1] */
    6 /* la matrice f ressort les dy */
    76
    87#include "titan.h"
    98
    10 void solve( double **a, double *f, int n0, int n1 )
     9void solve( double ***aa, int m, int n0, int n1 )
    1110{
    12    int    i,ii,imax,k,l,ll;
    13    double aamax,dum,*indx,sum,*vv,tmp;
     11   int    i,ii,imax,j,k,l,ll;
     12   double **a,aamax,dum,*indx,sum,**vv,tmp;
    1413   FILE   *out;
    1514
    1615   indx = dm1d( n0, n1 );
    17    vv   = dm1d( n0, n1 );
     16   vv   = dm2d( n0, n1, n0, n1 );
     17   a    = dm2d( n0, n1, n0, n1 );
    1818   imax = n0;
    1919
     20   for( i = n0; i <= n1; i++ ) for( j = n0; j <= n1; j++ ) a[i][j]=aa[m][i][j];
     21   
    2022   for( i = n0; i <= n1; i++ )
    2123   {
    … …  
    3234/*         aamax = 1.e-30; */
    3335      }
    34       vv[i] = 1.0e0 / aamax;   /* Save the scaling */
     36      vv[i][1] = 1.0e0 / aamax;   /* Save the scaling */
    3537/*
    3638      if( (aamax > 1.0e100)||(aamax < 1.0e-100) )
    … …  
    5961            sum -= ( a[i][l] * a[l][k] );
    6062         a[i][k] = sum;
    61          dum        = vv[i] * fabs(sum);  /* Figure of merit for the pivot */
     63         dum        = vv[i][1] * fabs(sum); /* Figure of merit for the pivot */
    6264         if( dum >= aamax )          /* Is it better than the best so far ? */
    6365         {
    … …  
    7476            a[k][l]    = dum;
    7577         }
    76          vv[imax] = vv[k];    /* Also interchange the scale factor */
     78         vv[imax][1] = vv[k][1];    /* Also interchange the scale factor */
    7779      }
    7880      indx[k] = imax;
    … …  
    9496   }                                  /* Go back to the next column in the reduction */
    9597
     98   for( i = n0; i <= n1; i++ )
     99   {
     100      for( k = n0; k <= n1; k++ )
     101         vv[i][k] = 0.0e0;
     102      vv[i][i] = 1.0e0;
     103   }
     104
     105   for( l = n0; l <= n1; l++ )
     106   {
    96107      ii = n0-1;
    97108      for( i = n0; i <= n1; i++ )
    98109      {
    99110         ll    = indx[i];
    100          sum   = f[ll];
    101          f[ll] = f[i];
     111         sum   = vv[ll][l];
     112         vv[ll][l] = vv[i][l];
    102113         if( ii != (n0-1) )
    103114            for( k = ii; k < i; k++ )
    104                sum -= ( a[i][k] * f[k] );
     115               sum -= ( a[i][k] * vv[k][l] );
    105116         else if( sum != 0.0e0 ) ii = i;
    106          f[i] = sum;
     117         vv[i][l] = sum;
    107118      }
    108119      for( i = n1; i >= n0; i-- )
    109120      {
    110          sum = f[i];
     121         sum = vv[i][l];
    111122         if( i < n1 )
    112123            for( k = i+1; k <= n1; k++ )
    113                sum -= ( a[i][k] * f[k] );
    114          f[i] = sum / a[i][i];
     124               sum -= ( a[i][k] * vv[k][l] );
     125         vv[i][l] = sum / a[i][i];
    115126      }
     127   }
    116128         
     129   for( i = n0; i <= n1; i++ )
     130     for( l = n0; l <= n1; l++ )
     131       aa[m][i][l]=vv[i][l];
     132
    117133   fdm1d( indx, n0 );
    118    fdm1d(   vv, n0 );
     134   fdm2d(   vv, n0, n1, n0 );
     135   fdm2d(    a, n0, n1, n0 );
    119136}

  • trunk/LMDZ.TITAN/libf/chimtitan/titan.h

    r3 r1126  
    77
    88#define R0    (double)(2575.0) /* Titan's radius */
    9 #define NLEV  (int)(55)   /* Nbre de niv verticaux - aussi dans titan_for.h */
     9#define NLEV  (int)(125)  /* Nbre de niv verticaux - =llm+70 dans common_mod */
     10#define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff */
     11#define NLRT  (int)(650)  /* Nbre de niv verticaux dans table fmoy - aussi dans common_mod */
     12
     13/* fluxes at 1300 km : upward is +, downward is - */
     14#define top_H   (double)(+1.1e4)
     15#define top_H2  (double)(+3.7e3)
     16#define top_N4S (double)(-1.1e8) /* = -2.5e8/2.27 ... */
    1017
    1118/* DEPEND DE LA VERSION CHIMIE: */
    1219#define VERCHIM "chimie_simpnit_051006_bis"
    13 #define NREAC (int)(377)    /* nombre de reactions - aussi dans titan_for.h */
    14 #define RDISS (int)(54)     /* nombre de photodiss - aussi dans titan_for.h */
    15 #define NC    (int)(44)     /* nb de composes      - aussi dans titan_for.h */
     20#define NREAC (int)(377)    /* nombre de reactions - aussi dans common_mod */
     21#define RDISS (int)(54)     /* nombre de photodiss - aussi dans common_mod */
     22#define NC    (int)(44)     /* nb de composes      - aussi dans common_mod */
    1623#define ST    (int)(NC)     /* nb de composes inverses */
    1724#define NHC   (int)(32)     /* nb hydrocarbons */
    … …  
    2734#endif
    2835
    29 /* void  gptitan_(const int *, float *, float *, float *, float *, float *,
    30       char (*)[10], float (*)[NLEV], float *,
    31       float *, float *, float *, int *, float *,
    32       float *, float (*)[NLEV][RDISS+1][15], float (*)[NLEV],
    33       int (*)[5], int *, int *, int (*)[200], int (*)[200][2] ); */
    3436void  chimie_(char (*)[10], double *, double *, double (*)[NLEV],
    3537              int (*)[5], int *, int *, int (*)[200][2], int (*)[200]);
    36 void  comp_(char (*)[10], double *);
    37 void  disso_(double (*)[NLEV][RDISS+1][15], int *);
    38 void  solve( double **, double *, int, int );
     38void  comp_(char (*)[10], double *, double *, double *, double (*)[NLEV]);
     39void  disso_(double (*)[NLRT][RDISS+1][15], int *);
     40double omega( double, double, double );
     41void  solve( double ***, int, int, int );
     42void  solve_b( double ***, double **, int, int, int );
    3943float *rm1d( int, int );
    4044float **rm2d( int, int, int, int );
Note: See TracChangeset for help on using the changeset viewer.