Changeset 1126 for trunk/LMDZ.TITAN/libf


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
Files:
3 added
2 deleted
14 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 );

  • trunk/LMDZ.TITAN/libf/phytitan/calchim.F

    r1056 r1126  
    11      SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
    2      .                   ctemp,cplay,cplev,
     2     .                   ctemp,cplay,cplev,czlay,czlev,
    33     .                   dqyc)
    44     
    … …  
    1010c            adaptation pour Titan 3D: 02/2009
    1111c            adaptation pour // : 04/2013
    12 c
     12c            extension chimie jusqu a 1300km : 10/2013
     13c
    1314c-------------------------------------------------
    1415c
    1516      use dimphy
    1617      use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
    17      .                     NLEV,NC,ND,NR
     18     .                     NLEV,NLRT,NC,ND,NR
    1819      USE comgeomphy,  only: rlatd
    19       use moyzon_mod,  only: klat
     20      USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
    2021      implicit none
    2122#include "dimensions.h"
    … …  
    3536      REAL         ctemp(nlon,klev)      ! Temperature
    3637      REAL         cplay(nlon,klev)      ! pression (Pa)
    37       REAL         cplev(nlon,klev)      ! pression intercouches (Pa)
     38      REAL         cplev(nlon,klev+1)    ! pression intercouches (Pa)
     39      REAL         czlay(nlon,klev)      ! altitude (m)
     40      REAL         czlev(nlon,klev+1)    ! altitude intercouches (m)
    3841     
    3942      REAL         dqyc(nlon,klev,NC)    ! Tendances especes chimiques
    … …  
    4649c variables envoyees dans la chimie: double precision
    4750
    48       REAL  temp_c(klev),press_c(klev)     ! T,p(mbar) a 1 lat donnee
    49       REAL  declin_c                       ! declinaison en degres
    50       REAL  cqy(klev,NC)                   ! frac mol qui seront modifiees
    51       REAL  surfhaze(klev)
    52       REAL  cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4)
    53       REAL  rinter(klev),nb(klev)
     51      REAL  temp_c(NLEV)
     52      REAL  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
     53      REAL  cqy(NLEV,NC)    ! frac mol qui seront modifiees
     54      REAL  cqy0(NLEV,NC)    ! frac mol avant modif
     55      REAL  surfhaze(NLEV)
     56      REAL  cprodaer(NLEV,4),cmaer(NLEV,4)
     57      REAL  ccsn(NLEV,4),ccsh(NLEV,4)
     58! rmil: milieu de couche, grille pour K,D,p,ct,T,y
     59! rinter: intercouche (grille RA dans la chimie)
     60      REAL  rmil(NLEV),rinter(NLEV),nb(NLEV)
     61      REAL,save :: kedd(NLEV)
     62
    5463      REAL  a,b
    5564      character str1*1,str2*2
    56       integer ntotftop
    57       parameter (ntotftop=50)
    58       integer nftop(ntotftop),isaisonflux
    59       REAL  fluxtop(NC),latit,tabletmp(ntotftop),factflux
    60       character*10 nomsftop(ntotftop+1)
     65      REAL  latit
    6166      character*20 formt1,formt2,formt3
    6267     
    … …  
    6772      SAVE    firstcal
    6873     
    69       REAL  mass(NC),duree
    70       REAL  tablefluxtop(NC,jjp1,5)
    71       REAL  botCH4
     74      integer dec,declinint,ialt
     75      REAL  declin_c                       ! declinaison en degres
     76      real  factalt,factdec,krpddec,krpddecp1,krpddecm1
     77      real  duree
     78      REAL,save :: mass(NC)
     79      REAL,save,allocatable :: md(:,:)
     80      REAL,save :: botCH4
    7281      DATA  botCH4/0.05/
     82      real,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
    7383      REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
    74       integer reactif(5,NR),nom_prod(NC),nom_perte(NC)
    75       integer prod(200,NC),perte(2,200,NC)
    76       SAVE    mass,tablefluxtop,botCH4
    77       SAVE    reactif,nom_prod,nom_perte,prod,perte
    78      
     84      integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
     85      integer,save :: prod(200,NC),perte(2,200,NC)
     86      character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
     87      real  ficalt,oldq,newq,zalt
     88      logical okfic
     89
    7990c-----------------------------------------------------------------------
    8091c***********************************************************************
    … …  
    94105c ************************************
    95106
    96         allocate(krpd(15,ND+1,klev,jjp1),krate(klev,NR))
    97 
    98 c Verification dimension verticale: coherence titan_for.h et klev
    99 c --------------------------------
    100 
    101         if (klev.ne.NLEV) then
    102            print*,'PROBLEME de coherence dans la dimension verticale:'
    103      .           ,klev,NLEV
    104            STOP
    105         endif
    106 
    107 c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro
     107        allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
     108
     109c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
    108110c ----------------------------------
    109111
    … …  
    114116        endif
    115117
    116 c calcul de temp_c, densites et press_c au milieu de l'ensemble des points:
    117 c ----------------------------------------------------------------------
    118 
    119         print*,'pression, densites et temp (chimie):'
     118c calcul de temp_c, densites et press_c en moyenne planetaire :
     119c --------------------------------------------------------------
     120
     121        print*,'pression, densites et temp (init chimie):'
    120122        print*,'level, press_c, nb, temp_c'
    121123        DO l=1,klev
     124         rinter(l)  = (zlevmoy(l)+RA)/1000.
     125         rmil(l)    = (zlaymoy(l)+RA)/1000.
    122126c     temp_c (K):
    123          temp_c(l)  = ctemp(nlon/2+1,l)
     127         temp_c(l)  = tmoy(l)
    124128c     press_c (mbar):
    125          press_c(l) = cplay(nlon/2+1,l)/100.
     129         press_c(l) = playmoy(l)/100.
    126130c     nb (cm-3):
    127131         nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
    128          print*,l,press_c(l),nb(l),temp_c(l)
     132         print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
    129133        ENDDO
     134
     135c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
     136       do l=klev+1,NLEV
     137         rinter(l) = rinter(klev)
     138     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     139         rmil(l)   = rmil(klev)
     140     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     141       enddo       
     142
     143c lecture de tcp.ver, une seule fois
     144c remplissage r1d,t1d,ct1d,p1d
     145        open(11,file='../../INPUT/tcp.ver',status='old')
     146        read(11,*) dummy
     147        do i=1,131
     148          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
     149c         print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
     150        enddo
     151        close(11)
     152
     153c extension pour klev+1 a NLEV avec tcp.ver
     154c la jonction klev/klev+1 est brutale... Tant pis ?
     155        ialt=1
     156        do l=klev+1,NLEV
     157           zalt=rmil(l)-RA/1000.
     158           do i=ialt,130
     159            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
     160              ialt=i
     161            endif
     162           enddo
     163           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
     164           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)
     165     &                      + log(p1d(ialt+1)) * factalt    )
     166           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt)
     167     &                      + log(ct1d(ialt+1)) * factalt    )
     168           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
     169           print*,l,zalt,press_c(l),nb(l),temp_c(l)
     170        enddo
    130171       
    131172c caracteristiques des composes:       
    132173c -----------------------------
    133         mass = 0.0
    134         call comp(nomqy_c,mass)
     174        mass(:) = 0.0
     175        call comp(nomqy_c,nb,temp_c,mass,md)
    135176        print*,'           Mass'
    136177        do ic=1,NC
    137178          print*,nomqy_c(ic),mass(ic)
     179c         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
    138180        enddo
    139181       
    140182
    141 c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON...
    142 
    143 c flux dans la derniere couche:(issu du modele 1D, a 500 km)
    144 c -----------------------------
    145 c       !!  en cm-2.s-1 !!
    146 c  ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54
    147 c
    148 c Lecture des tables de flux en fonction lat et saison
    149 c ----------------------------------------------------
    150 
    151         WRITE(str2(1:2),'(i2.2)') ntotftop
    152         formt1 = '(A10,'//str2//'(3X,A10))'
    153         formt2 = '(F12.6,'//str2//'(2X,F13.6))'
    154         formt3 = '(F13.6,'//str2//'(2X,F13.6))'
    155        
    156         do j=1,jjp1
    157           do ic=1,NC
    158             do l=1,5
    159              tablefluxtop(ic,j,l) = 0.
    160             enddo
    161           enddo
    162         enddo
    163        
    164         do l=1,4
    165           WRITE(str1(1:1),'(i1.1)') l
    166 c hx1 -> Ls=RPI/2 / hx4 -> Ls=0
    167           open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1)
    168 c         open(10,file="FLUXTOP/flux500.hs"//str1)
    169 c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ...
    170 
    171           read(10,formt1) nomsftop
    172           do j=1,ntotftop
    173            do i=1,10   !justification a gauche...
    174              if (nomsftop(j+1)(i:i).ne.' ') then
    175                 nomsftop(j) = nomsftop(j+1)(i:)
    176                 goto 100
    177              endif
    178            enddo
    179 100        continue
    180 c         print*,nomsftop(j)
    181            nftop(j) = 0
    182            do i=1,NC
    183             if (nomqy_c(i).eq.nomsftop(j)) then
    184              nftop(j) = i
    185             endif
    186            enddo
    187 c          if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j))
    188           enddo
    189 
    190          
    191 c lecture des flux. Les formats sont un peu alambiques...
    192 c     a revoir eventuellement
    193           do j=1,jjp1/2+1
    194             read(10,formt2) latit,tabletmp
    195             do i=1,ntotftop
    196               if (nftop(i).ne.0) then
    197                tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
    198               endif
    199             enddo
    200           enddo
    201           do j=jjp1/2+2,jjp1
    202             read(10,formt3) latit,tabletmp
    203             do i=1,ntotftop
    204               if (nftop(i).ne.0) then
    205                tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
    206               endif
    207             enddo
    208           enddo
    209 
    210           close(10)
    211         enddo  ! l
    212        
    213         do j=1,jjp1
    214           do ic=1,NC
    215              tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5)
    216           enddo
    217         enddo
    218        
    219 c  Stockage des composes utilises dans la prod d'aerosols
     183c  Stockage des composes utilises dans la prod d aerosols
    220184c     (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
    221185c     ! decalage de 1 car utilise dans le c !
    222 c ------------------------------------------
    223 c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN
    224 c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10
    225 c ------------------------------------------
    226186
    227187        do ic=1,NC
    … …  
    244204          if (nomqy_c(ic).eq."HCN") then
    245205            utilaer(6)  = ic-1
    246             do j=1,jjp1
    247               do i=1,5
    248              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6.
    249               enddo
    250             enddo
    251206          endif
    252207          if (nomqy_c(ic).eq."HC3N") then
    … …  
    268223          if (nomqy_c(ic).eq."C2H2") then
    269224            utilaer(3)  = ic-1
    270             do j=1,jjp1
    271               do i=1,5
    272 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3
    273              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10.
    274               enddo
    275             enddo
    276225          endif
    277226          if (nomqy_c(ic).eq."AC6H6") then
    … …  
    298247          if (nomqy_c(ic).eq."AC6H5") then
    299248            utilaer(12) = ic-1
    300           endif
    301 
    302           if (nomqy_c(ic).eq."C2H6") then
    303             do j=1,jjp1
    304               do i=1,5
    305 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640.
    306              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200.
    307               enddo
    308             enddo
    309           endif
    310           if (nomqy_c(ic).eq."C3H8") then
    311             do j=1,jjp1
    312               do i=1,5
    313 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
    314              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.)
    315               enddo
    316             enddo
    317           endif
    318           if (nomqy_c(ic).eq."C4H10") then
    319             do j=1,jjp1
    320               do i=1,5
    321              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
    322               enddo
    323             enddo
    324249          endif
    325250
    … …  
    369294c       enddo
    370295
     296
     297c   init kedd
     298c   ---------
     299c   kedd stays constant with time and space
     300c   => linked to pressure rather than altitude...
     301 
     302      kedd(:) = 5e2 ! valeur mise par defaut 
     303               ! UNITE ? doit etre ok pour gptitan
     304      do l=1,NLEV
     305       zalt=rmil(l)-RA/1000.  ! z en km
     306       if     ((zalt.ge.250.).and.(zalt.le.400.)) then
     307         kedd(l) = 10.**(3.+(zalt-250.)/50.)
     308! 1E3 at 250 km
     309! 1E6 at 400 km
     310       elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
     311         kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.)
     312! 2E7 at 600 km
     313       elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
     314         kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.)
     315! 1E8 above 900 km
     316       elseif ( zalt.gt.900.                    ) then
     317        kedd(l) = 1.e8
     318       endif
     319      enddo
     320
    371321      ENDIF  ! premier appel
    372322
    … …  
    380330c      print*,'declinaison en degre=',declin_c
    381331       
    382 c-----------------------------------------------------------------------
    383 c
    384 c   calcul du facteur d'interpolation entre deux saisons pour flux
    385 c   --------------------------------------------------------------
    386 
    387        isaisonflux = int(ls_rad*2./RPI)+1
    388        if (isaisonflux.eq.5) then
    389          isaisonflux = 1
    390        endif
    391        factflux =    (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/
    392      .   (sin(isaisonflux*RPI/2.)-sin((isaisonflux-1)*RPI/2.))
    393 
    394332c***********************************************************************
    395333c***********************************************************************
    … …  
    410348c***********************************************************************
    411349
     350c-----------------------------------------------------------------------
     351c
     352c   Distance radiale (intercouches, en km)
     353c   ----------------------------------------
     354
     355       do l=1,klev
     356         rinter(l) = (RA+czlev(j,l))/1000.
     357         rmil(l)   = (RA+czlay(j,l))/1000.
     358c        print*,rinter(l)
     359       enddo
     360
     361c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
     362       do l=klev+1,NLEV
     363         rinter(l) = rinter(klev)
     364     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     365         rmil(l)   = rmil(klev)
     366     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     367       enddo
     368       
    412369c-----------------------------------------------------------------------
    413370c
    … …  
    423380               nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
    424381       ENDDO
    425        
    426 c-----------------------------------------------------------------------
    427 c
    428 c   Distances radiales (intercouches, en km)
    429 c   ----------------------------------------
    430 
    431        rinter(1) = RA/1000. 
    432        do l=2,klev
    433 c A REVOIR: g doit varier avec r !
    434          rinter(l) = rinter(l-1) +
    435      .    (cplev(j,l-1)-cplev(j,l))/cplay(j,l)*(RD*temp_c(l)/RG)/1000.
    436 c        print*,rinter(l)
    437        enddo
    438 
    439 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    440 c FLUX AU TOP: interpolation en saison et
    441 c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C)
    442 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    443         do ic=1,NC
    444           fluxtop(ic) = tablefluxtop(ic,j,isaisonflux)  *(1-factflux)
    445      .                + tablefluxtop(ic,j,isaisonflux+1)*  factflux
    446           fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5)
    447         enddo
    448 
    449 c TEST: sortie de l'un des flux avec saveLs et factflux
    450 c       dans un fichier special pour tracer avec gnuplot
    451 c       if (j.eq.2) then
    452 c       open(unit=11,file="flux_80N.txt",status='old',position='append')
    453 c         write(11,*) ls_rad*180./RPI, factflux,
    454 c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
    455 c         write(11,*) ls_rad*180./RPI, factflux,
    456 c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
    457 c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
    458 c       close(11)
    459 c       endif
    460 c       if (j.eq.17) then
    461 c       open(unit=11,file="flux_eq.txt",status='old',position='append')
    462 c         write(11,*) ls_rad*180./RPI, factflux,
    463 c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
    464 c         write(11,*) ls_rad*180./RPI, factflux,
    465 c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
    466 c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
    467 c       close(11)
    468 c       endif
    469 c FIN TEST: sortie de l'un des flux avec ls et factflux
    470 
    471 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     382c extension pour klev+1 a NLEV avec tcp.ver
     383       ialt=1
     384       do l=klev+1,NLEV
     385           zalt=rmil(l)-RA/1000.
     386           do i=ialt,130
     387            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
     388              ialt=i
     389            endif
     390           enddo
     391           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
     392           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)
     393     &                      + log(p1d(ialt+1)) * factalt    )
     394           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt)
     395     &                      + log(ct1d(ialt+1)) * factalt    )
     396           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
     397       enddo
     398               
     399c-----------------------------------------------------------------------
     400c
     401c   passage krpd => krate
     402c   ---------------------
     403c initialisation krate pour dissociations
     404
     405      if ((declin_c*10+267).lt.14.) then
     406          declinint = 0
     407          dec       = 0
     408      else
     409       if ((declin_c*10+267).gt.520.) then
     410          declinint = 14
     411          dec       = 534
     412       else
     413          declinint = 1
     414          dec       = 27
     415          do while( (declin_c*10+267).ge.real(dec+20) )
     416            dec       = dec+40
     417            declinint = declinint+1
     418          enddo
     419       endif
     420      endif
     421      if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
     422          factdec = ( declin_c - (dec-267)/10. ) / 4.
     423      else
     424          factdec = ( declin_c - (dec-267)/10. ) / 2.7
     425      endif
     426
     427      do l=1,NLEV
     428
     429c INTERPOL EN ALT POUR k (krpd tous les deux km)
     430        ialt    = int((rmil(l)-RA/1000.)/2.)+1
     431        factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
     432
     433        do i=1,ND+1
     434          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt)
     435     &              + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
     436          krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt)
     437     &              + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
     438          krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt)
     439     &              + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
     440
     441                    ! ND+1 correspond a la dissociation de N2 par les GCR
     442          if ( factdec.lt.0. ) then
     443        krate(l,i) = krpddecm1 * abs(factdec)
     444     &             + krpddec   * ( 1 + factdec)
     445          endif
     446          if ( factdec.gt.0. ) then
     447        krate(l,i) = krpddecp1 * factdec
     448     &             + krpddec   * ( 1 - factdec)
     449          endif
     450          if ( factdec.eq.0. ) krate(l,i) = krpddec
     451        enddo       
     452      enddo       
    472453
    473454c-----------------------------------------------------------------------
    … …  
    475456c   composition
    476457c   ------------
    477        
     458
    478459       do ic=1,NC
    479460        do l=1,klev
    480           cqy(l,ic) = qy_c(j,l,ic)
    481         enddo
    482        enddo
    483              
     461          cqy(l,ic)  = qy_c(j,l,ic)
     462          cqy0(l,ic) = cqy(l,ic)
     463        enddo
     464       enddo
     465
     466c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
     467
     468      WRITE(str2,'(i2.2)') klat(j)
     469      fich = "comp_"//str2//".dat"
     470      inquire (file=fich,exist=okfic)
     471      if (okfic) then
     472       open(11,file=fich,status='old')
     473c premiere ligne=declin
     474       read(11,'(A15)') ficdec
     475       write(curdec,'(E15.5)') declin_c
     476c si la declin est la meme: ce fichier a deja ete reecrit
     477c on lit la colonne t-1 au lieu de la colonne t
     478c (cas d une bande de latitude partagee par 2 procs)
     479       do ic=1,NC
     480        read(11,'(A10)') name
     481        if (name.ne.nomqy_c(ic)) then
     482          print*,"probleme lecture ",fich
     483          print*,name,nomqy_c(ic)
     484          stop
     485        endif
     486        if (ficdec.eq.curdec) then
     487          do l=klev+1,NLEV
     488            read(11,*) ficalt,cqy(l,ic),newq
     489          enddo
     490        else
     491          do l=klev+1,NLEV
     492            read(11,*) ficalt,oldq,cqy(l,ic)
     493          enddo
     494        endif
     495       enddo
     496       close(11)
     497      else       ! le fichier n'est pas la
     498       do ic=1,NC
     499        do l=klev+1,NLEV
     500          cqy(l,ic)=cqy(klev,ic)
     501        enddo
     502       enddo
     503      endif
     504      cqy0 = cqy
     505
    484506c-----------------------------------------------------------------------
    485507c
    … …  
    502524       
    503525       call gptitan(rinter,temp_c,nb,
    504      $              nomqy_c,cqy,fluxtop,
    505      $              declin_c,duree,(klat(j)-1),mass,
    506      $              botCH4,krpd,krate,reactif,
     526     $              nomqy_c,cqy,
     527     $              duree,(klat(j)-1),mass,md,
     528     $              kedd,botCH4,krate,reactif,
    507529     $              nom_prod,nom_perte,prod,perte,
    508530     $              aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
    … …  
    514536       do ic=1,NC
    515537        do l=1,klev
    516           dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim  ! en /s
     538          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
    517539        enddo
    518540       enddo
    … …  
    535557
    536558       endif
     559
     560c-----------------------------------------------------------------------
     561c
     562c   sauvegarde compo sur NLEV
     563c   -------------------------
     564
     565c dans fichier compo_klat(j) (01 à 49)
     566       
     567      open(11,file=fich,status='replace')
     568c premiere ligne=declin
     569      write(11,'(E15.5)') declin_c
     570      do ic=1,NC
     571       write(11,'(A10)') nomqy_c(ic)
     572       do l=klev+1,NLEV
     573        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
     574     .                                cqy0(l,ic),cqy(l,ic)
     575       enddo
     576      enddo
     577      close(11)
    537578
    538579c***********************************************************************

  • trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90

    r1056 r1126  
    3434
    3535! ancien titan_for.h
    36       INTEGER, PARAMETER :: NLEV=llm,NC=44,ND=54,NR=377
     36      INTEGER, PARAMETER :: NLEV=llm+70,NC=44,ND=54,NR=377
    3737!$OMP THREADPRIVATE(NLEV,NC,ND,NR)
    3838!!!  doivent etre en accord avec titan.h
     39! pour l'UV (650 niveaux de 2 km)
     40      INTEGER, PARAMETER :: NLRT=650
    3941
    4042! ancien diagmuphy.h

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F

    r1080 r1126  
    162162          XICLDI2(K)=TNI
    163163 420  CONTINUE
     164
     165c DEBUG
     166c       print*,"wnoi=",WNOI
    164167
    165168C
    … …  
    284287        TauGID(ig,:,:) = TAUGID_1pt(:,:)
    285288
     289c DEBUG
     290c     if(ig.eq.(ngrid/2+16)) then
     291c         print*,ig,'/',KLON,':'
     292c         print*,'TauHID_1',TAUHID(ig,1,:)
     293c         print*,'TauGID_1',TAUGID(ig,1,:)
     294c         print*,'TauHID_50',TAUHID(ig,50,:)
     295c         print*,'TauGID_50',TAUGID(ig,50,:)
     296c         print*,"DTAUI_1=",DTAUI(ig,1,:)
     297c         print*,"DTAUI_50=",DTAUI(ig,50,:)
     298c         print*,'cosBI_1',COSBI(ig,1,:)
     299c         print*,'cosBI_50',COSBI(ig,50,:)
     300c         print*,'WBARI_1',WBARI(ig,1,:)
     301c         print*,'WBARI_50',WBARI(ig,50,:)
     302c         stop
     303c     endif
     304
    286305c************************************************************************
    287306c************************************************************************

  • trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F

    r1058 r1126  
    282282c     CBAR=xv3(j,k)
    283283
     284c     print*, 'HERE, CIRS AEROSOLS'
     285c     call cirs_haze(PRESS(J),WNOI(K),TAEROS,TAEROSCAT,CBAR)
    284286
    285287c*********** EN TRAVAUX ***************************

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F

    r1081 r1126  
    148148c       close(1)
    149149
     150c DEBUG
     151c       print*,"wnov=",WNOV
     152
    150153      endif    ! fin initialisations premier appel
    151154
    … …  
    167170        else
    168171         if (ig.eq.1)  then
    169 c initialisation zqaer_1pt a partir d'une look-up table (uniforme en ig)
     172c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
    170173c boucle sur nrad=10
    171174           open(10,file="qaer_eq_1d.dat")
    … …  
    220223        TauGVD(ig,:,:) = TAUGVD_1pt(:,:)
    221224
     225c DEBUG
     226c     if(ig.eq.(ngrid/2+16)) then
     227c         print*,ig,'/',KLON,':'
     228c         print*,'TauHVD_1',TAUHVD(ig,1,:)
     229c         print*,'TauGVD_1',TAUGVD(ig,1,:)
     230c         print*,'TauHVD_50',TAUHVD(ig,50,:)
     231c         print*,'TauGVD_50',TAUGVD(ig,50,:)
     232c     stop
     233c     endif
     234
    222235 101  CONTINUE
    223236

  • trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F

    r1058 r1126  
    231231
    232232       IF (TAEROSCAT.le.0.) CBAR=0.
     233
     234c     print*, 'HERE, CIRS AEROSOLS'
     235c     call cirs_haze(PRESS(J),WNOV(K),TAEROS,TAEROSCAT,CBAR)
    233236
    2342371699  FORMAT(a3,2I3,3(ES15.7,1X))

  • trunk/LMDZ.TITAN/libf/phytitan/physiq.F

    r1120 r1126  
    419419         itap    = 0
    420420         itaprad = 0
    421          itapchim    = 0
     421         itapchim    = 1
    422422
    423423c init rnuabar
    … …  
    11261126         tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
    11271127     .                        + d_tr_mph(:,:,1:nmicro)*dtime
    1128 
     1128c        call WriteField_phy('physiq_d_tr_mph01',
     1129c    .                          d_tr_mph(:,:,1),klev)
     1130c        call WriteField_phy('physiq_d_tr_mph10',
     1131c    .                          d_tr_mph(:,:,10),klev)
     1132        endif
    11291133c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
    11301134c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
    … …  
    11411145        ENDDO
    11421146
    1143         endif
    1144 
    11451147c condensation:
    11461148c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
    … …  
    11491151     .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
    11501152        endif
     1153
     1154c chimie:
    11511155        if ((chimi).and.(nqmax.gt.nmicro)) then
    1152 c chimie:
    11531156         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
    11541157        endif
    11551158
    1156       endif
     1159      endif  !iflag_trac=1
    11571160
    11581161c       ch4=0.

  • trunk/LMDZ.TITAN/libf/phytitan/phytrac.F

    r1072 r1126  
    459459c ------------
    460460       CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
    461      .              ztemp,zplay,zplev,
     461     .              ztemp,zplay,zplev,zzlay,zzlev,
    462462     .              pdyfi)   
    463463c ychim ne doit pas etre modifie, pdyfi en /s

  • trunk/LMDZ.TITAN/libf/phytitan/profile.F

    r904 r1126  
    9494       b2 =       3183.
    9595       c2 =       4737.
     96c pres est en Pa => conversion car expression veut p en hPa
    9697       DO il=1,nlev
    97          temp(il)=a1*exp(-((pres(il)-b1)/c1)**2.)
    98      .          + a2*exp(-((pres(il)-b2)/c2)**2.)
     98         temp(il)=a1*exp(-((pres(il)/100.-b1)/c1)**2.)
     99     .          + a2*exp(-((pres(il)/100.-b2)/c2)**2.)
    99100       ENDDO
    100101       zkm(1)  = 0.0
Note: See TracChangeset for help on using the changeset viewer.