Changeset 1956


Ignore:
Timestamp:
Jul 1, 2018, 5:07:48 PM (7 years ago)
Author:
jvatant
Message:

+ Add surface methane CH4 flux pseudo-evap diagnostic and include it in bottom layer zdqcond as >0.
+ minor cosmetic changes in gptitan.c
--JVO

Location:
trunk/LMDZ.TITAN
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/README

    r1950 r1956  
    14791479Optimization for chemistry firstcall C routine disso.c
    14801480-> Row-major-friendly declaration and nesting of loops for huge array krpd
    1481 -> ~factor 4 for this routine who was ~ 1/2 of firstcall -> ~ gain ~5% total run-time for standard runs
     1481-> ~factor 4 for this routine who was ~ 1/2 of firstcall -> ~ gain ~5% total run-time for standard runs
     1482
     1483== 01/07/18 == JVO
     1484Add surface methane CH4 flux pseudo-evap diagnostic and include it in bottom layer zdqcond as >0.

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

    r1908 r1956  
    99
    1010void gptitan_(
    11      double *RA, double *TEMP, double *NB,
    12      char CORPS[][10], double Y[][NLEV],
    13      double *FIN, double *LAT, double *MASS, double MD[][NLEV],
    14      double *KEDD, double *botCH4, double KRATE[][NLEV],
    15      int reactif[][5], int *nom_prod, int *nom_perte,
    16      int prod[][200], int perte[][200][2], int *aerprod, int *utilaer,
    17      double MAER[][NLEV], double PRODAER[][NLEV],
    18      double CSN[][NLEV], double CSH[][NLEV],
    19      int *htoh2, double *surfhaze)
     11              double *RA, double *TEMP, double *NB,
     12              char CORPS[][10], double Y[][NLEV],
     13              double *FIN, double *LAT, double *MASS, double MD[][NLEV],
     14              double *KEDD, double *botCH4, double *fluxCH4, double KRATE[][NLEV],
     15              int reactif[][5], int *nom_prod, int *nom_perte,
     16              int prod[][200], int perte[][200][2], int *aerprod, int *utilaer,
     17              double MAER[][NLEV], double PRODAER[][NLEV],
     18              double CSN[][NLEV], double CSH[][NLEV],
     19              int *htoh2, double *surfhaze)
    2020{
    21    char   outlog[100],corps[100][10];
    22    int    i,j,k,l;
    23    int    ireac,ncom1,ncom2;
    24    double  ***a,***b,**c;
    25    double  *fl,*fp,*mu,**jac,**ym1,**f;
    26    double  fluxCH4;
    27    double  conv,delta,deltamax;
    28    double  cm,cp,dim,dip,dm,dp,dym,dyp,km,kp,r,dra,dram,drap;
    29    double  np,nm,s,test,time,ts,v,dv;
    30    char   str2[15];
    31    FILE   *out;
    32 
    33 /* va avec htoh2 */
    34    double  dyh,dyh2;
    35 
    36 /* va avec aer */
    37    double  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
    38    double  **k_dep,**faer;
    39    double  *productaer,*csurn,*csurh,*mmolaer;
    40 
    41    if( (*aerprod) == 1 )
    42    {
    43     k_dep = dm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
    44     faer  = dm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
    45     productaer  = dm1d( 0, 3 );     /* local production rate by pathways */
    46     mmolaer     = dm1d( 0, 3 );     /* local molar mass by pathways */
    47     csurn = dm1d( 0, 3 );           /* local C/N by pathways */
    48     csurh = dm1d( 0, 3 );           /* local C/H by pathways */
    49    }
    50 
    51 /* DEBUG */
    52       printf("CHIMIE: lat=%g\n",(*LAT));
    53 /**/
    54 
    55    for( i = 0; i <= NC; i++)
    56    {
    57      strcpy( corps[i], CORPS[i] );
    58      corps[i][strcspn(CORPS[i], " ")] = '\0';
    59    }
     21  char   outlog[100],corps[100][10];
     22  int    i,j,k,l;
     23  int    ireac,ncom1,ncom2;
     24  int    in2, ich4, ih, ih2, in4s;
     25  double  ***a,***b,**c;
     26  double  *fl,*fp,*mu,**jac,**ym1,**f;
     27  double  dyCH4surf;
     28  double  conv,delta,deltamax;
     29  double  cm,cp,dim,dip,dm,dp,dym,dyp,km,kp,r,dra,dram,drap;
     30  double  np,nm,s,test,time,ts,v,dv;
     31  char   str2[15];
     32  FILE   *out;
     33
     34  /* va avec htoh2 */
     35  double  dyh,dyh2;
     36
     37  /* va avec aer */
     38  double  dyc2h2,dyhc3n,dyhcn,dynccn,dych3cn,dyc2h3cn;
     39  double  **k_dep,**faer;
     40  double  *productaer,*csurn,*csurh,*mmolaer;
     41
     42  if( (*aerprod) == 1 )
     43    {
     44      k_dep = dm2d( 1, 5, 1, 3 );     /* k en s-1, reactions d'initiation */
     45      faer  = dm2d( 1, 5, 1, 3 );     /* fraction de chaque compose */
     46      productaer  = dm1d( 0, 3 );     /* local production rate by pathways */
     47      mmolaer     = dm1d( 0, 3 );     /* local molar mass by pathways */
     48      csurn = dm1d( 0, 3 );           /* local C/N by pathways */
     49      csurh = dm1d( 0, 3 );           /* local C/H by pathways */
     50    }
     51
     52  /* DEBUG */
     53  printf("CHIMIE: lat=%g\n",(*LAT));
     54  /**/
     55
     56  for( i = 0; i <= NC; i++)
     57    {
     58      strcpy( corps[i], CORPS[i] );
     59      corps[i][strcspn(CORPS[i], " ")] = '\0';
     60    }
    6061   
    61    strcpy( outlog, "chimietitan" );
    62    strcat( outlog, ".log" );
    63    out = fopen( outlog, "w" );
    64    fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
    65    fclose( out );
    66 
    67    deltamax = 1.e5;
    68    test = 1.0e-15;
    69 
    70 /* valeur de r:
    71 r = g0 R0^2 / R * 2 * 1E-3
    72 avec g0 en cm/s2, R0 en km, mu et mass en g
    73 */
    74    r        = 21.595656e0;
    75 
    76 /*  DEBUG
    77             out = fopen( outlog, "a" );
    78    fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
    79             fclose( out );
    80 */
    81    fl      = dm1d( 0,   NC-1 );
    82    fp      = dm1d( 0,   NC-1 );
    83    mu      = dm1d( 0, NLEV-1 );
    84    ym1     = dm2d( 0,   NC-1, 0, NLEV-1 );
    85    f       = dm2d( 0,   NC-1, 0, NLEV-1 );
    86    jac     = dm2d( 0,   NC-1, 0,   NC-1 );
    87    c       = dm2d( 0, NLEV-1, 0,   NC-1 );
    88    a       = dm3d( 0, NLEV-1, 0,   NC-1, 0,   NC-1 );
    89    b       = dm3d( 0, NLEV-1, 0,   NC-1, 1,    2 );
    90 
    91 /* DEBUG */
    92 /*
    93             out = fopen( "err.log", "a" );
    94             fprintf( out,"%s\n", );
    95             fclose( out );
    96 */
    97 
    98 /* initialisation mu, CH4 au sol */
     62  strcpy( outlog, "chimietitan" );
     63  strcat( outlog, ".log" );
     64  out = fopen( outlog, "w" );
     65  fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
     66  fclose( out );
     67
     68  deltamax = 1.e5;
     69  test = 1.0e-15;
     70
     71  /* valeur de r:
     72     r = g0 R0^2 / R * 2 * 1E-3
     73     avec g0 en cm/s2, R0 en km, mu et mass en g
     74  */
     75  r        = 21.595656e0;
     76
     77  /*  DEBUG
     78      out = fopen( outlog, "a" );
     79      fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
     80      fclose( out );
     81  */
     82  fl      = dm1d( 0,   NC-1 );
     83  fp      = dm1d( 0,   NC-1 );
     84  mu      = dm1d( 0, NLEV-1 );
     85  ym1     = dm2d( 0,   NC-1, 0, NLEV-1 );
     86  f       = dm2d( 0,   NC-1, 0, NLEV-1 );
     87  jac     = dm2d( 0,   NC-1, 0,   NC-1 );
     88  c       = dm2d( 0, NLEV-1, 0,   NC-1 );
     89  a       = dm3d( 0, NLEV-1, 0,   NC-1, 0,   NC-1 );
     90  b       = dm3d( 0, NLEV-1, 0,   NC-1, 1,    2 );
     91
     92  /* DEBUG */
     93  /*
     94     out = fopen( "err.log", "a" );
     95     fprintf( out,"%s\n", );
     96     fclose( out );
     97  */
     98
     99  /* init indexes of special species */
     100  for( i = 0; i <= ST-1; i++ )
     101    {
     102      if( strcmp(corps[i], "CH4") == 0 ) ich4 = i;
     103      if( strcmp(corps[i], "H"  ) == 0 ) ih   = i;
     104      if( strcmp(corps[i], "H2" ) == 0 ) ih2  = i;
     105      if( strcmp(corps[i], "N2" ) == 0 ) in2  = i;
     106      if( strcmp(corps[i], "N4S") == 0 ) in4s = i;
     107    }
     108
     109  /* initialisation mu, CH4 au sol */
    99110     
    100    for( j = 0; j <= NLEV-1; j++ )
    101    {
     111  for( j = 0; j <= NLEV-1; j++ )
     112    {
    102113      mu[j] = 0.0e0;
    103114      for( i = 0; i <= ST-1; i++ )
    104       {
    105          if( ( strcmp(corps[i], "CH4") == 0 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
    106          {
    107              fluxCH4 = (*botCH4 - Y[i][j]);
    108              Y[i][j] = *botCH4;
    109          }
    110          mu[j] += ( MASS[i] * Y[i][j] );
    111       }
    112    }
    113 
    114 /* initialisation compo avant calcul */
    115    for( j = NLEV-1; j >= 0; j-- )
    116       for( i = 0; i <= ST-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
    117 
    118 /*
    119 ==========================================================================
    120    STRATEGIE:
    121     INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
    122     PUIS INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
    123 ==========================================================================
    124 
    125 PREMIERE ETAPE:
    126 ===============
    127     INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
    128 ===============
    129 */
    130 
    131 /* ****************** */
    132 /*  Time loop:        */
    133 /* ****************** */
    134 
    135 time     = ts = 0.0e0;
    136 delta    = 1.e-3;
    137 
    138 while( time < (*FIN) )     
    139 {
    140 
    141 
    142 /* DEBUG
    143    for( j = NLEV-1; j >= NLD; j-- )
    144    {
    145             out = fopen( outlog, "a" );
    146         fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
    147             fclose( out );
    148 
    149             out = fopen( "profils.log", "a" );
    150     fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
    151     for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
    152     for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
    153             fclose( out );
     115        {
     116          if( ( i == ich4 ) && ( Y[i][j] <= *botCH4 ) && ( j == 0 ) )
     117            {
     118              dyCH4surf = (*botCH4 - Y[i][j]);
     119              Y[i][j] = *botCH4;
     120            }
     121          mu[j] += ( MASS[i] * Y[i][j] );
     122        }
    154123    }
    155    exit(0);
    156 */
    157 
    158 
    159 /* ------------------------------ */
    160 /* Calculs variations et jacobien */
    161 /* ------------------------------ */
    162 
    163    for( j = NLEV-1; j >= NLD; j-- )
    164    {
    165 
    166 /* init of step */
    167 /* ------------ */
    168          for( i = 0; i <= ST-1; i++ )
    169          {
    170             fp[i] = fl[i] = 0.0e0;
    171             for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
    172          }
    173 
    174 /* Chimie */
    175 /* ------ */
    176 
    177 /* productions et pertes chimiques */
    178          for( i = 0; i <= ST-1; i++ )   
    179          {
    180             Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
    181 
    182             for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
    183             {
    184                ireac = prod[i][l];                  /* Number of the reaction involves. */
    185                ncom1 = reactif[ireac][0];           /* First compound which reacts. */
    186                if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
    187                {
    188                   jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
    189                   fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
    190                }
    191                else                                 /* General case. */
    192                {
    193                   ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
    194                   jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
    195                   jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
    196                   fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
    197                }
    198             }
     124
     125  /* initialisation compo avant calcul */
     126  for( j = NLEV-1; j >= 0; j-- )
     127    for( i = 0; i <= ST-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
     128
     129  /*
     130     ==========================================================================
     131     STRATEGIE:
     132     INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
     133     PUIS INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
     134     ==========================================================================
     135
     136     PREMIERE ETAPE:
     137     ===============
     138     INVERSION COMPLETE AVEC DIFFUSION ENTRE NLEV-1 et NLD
     139     ===============
     140  */
     141
     142  /* ****************** */
     143  /*  Time loop:        */
     144  /* ****************** */
     145
     146  time     = ts = 0.0e0;
     147  delta    = 1.e-3;
     148
     149  while( time < (*FIN) )     
     150    {
     151
     152
     153      /* DEBUG
     154         for( j = NLEV-1; j >= NLD; j-- )
     155         {
     156         out = fopen( outlog, "a" );
     157         fprintf(out,"j=%d z=%e nb=%e T=%e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
     158         fclose( out );
     159
     160         out = fopen( "profils.log", "a" );
     161         fprintf(out,"%d %e %e %e\n",j,(RA[j]-R0),NB[j],TEMP[j]);
     162         for (i=0;i<=NREAC-1;i++) fprintf(out,"%d %e\n",i,KRATE[i][j]);
     163         for (i=0;i<=ST-1;i++) fprintf(out,"%10s %e\n",corps[i],Y[i][j]);
     164         fclose( out );
     165         }
     166         exit(0);
     167      */
     168
     169
     170      /* ------------------------------ */
     171      /* Calculs variations et jacobien */
     172      /* ------------------------------ */
     173
     174      for( j = NLEV-1; j >= NLD; j-- )
     175        {
     176
     177          /* init of step */
     178          /* ------------ */
     179          for( i = 0; i <= ST-1; i++ )
     180            {
     181              fp[i] = fl[i] = 0.0e0;
     182              for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
     183            }
     184
     185          /* Chimie */
     186          /* ------ */
     187
     188          /* productions et pertes chimiques */
     189          for( i = 0; i <= ST-1; i++ )   
     190            {
     191              Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
     192
     193              for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
     194                {
     195                  ireac = prod[i][l];                  /* Number of the reaction involves. */
     196                  ncom1 = reactif[ireac][0];           /* First compound which reacts. */
     197                  if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     198                    {
     199                      jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
     200                      fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
     201                    }
     202                  else                                 /* General case. */
     203                    {
     204                      ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
     205                      jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
     206                      jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
     207                      fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
     208                    }
     209                }
    199210           
    200             for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
    201             {
    202                ireac = perte[i][l][0];              /* Reaction number. */
    203                ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
    204                if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
    205                {
    206                   jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
    207                   fl[i]     += ( KRATE[ireac][j] * NB[j] );
    208                }
    209                else                                 /* General case. */
    210                {
    211                   jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
    212                   jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
    213                   fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
    214                }
    215             }
    216          }
    217 
    218 
    219 /* Aerosols */
    220 /* -------- */
    221          if( (*aerprod) == 1 )
    222          {
    223              aer(corps,TEMP,NB,Y,&j,k_dep,faer,
    224               &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
    225               mmolaer,productaer,csurn,csurh);
    226 
    227              for( i = 0; i <= 3; i++ )
    228              {
    229                PRODAER[i][j] = productaer[i];
     211              for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
     212                {
     213                  ireac = perte[i][l][0];              /* Reaction number. */
     214                  ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
     215                  if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     216                    {
     217                      jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
     218                      fl[i]     += ( KRATE[ireac][j] * NB[j] );
     219                    }
     220                  else                                 /* General case. */
     221                    {
     222                      jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
     223                      jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
     224                      fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
     225                    }
     226                }
     227            }
     228
     229
     230          /* Aerosols */
     231          /* -------- */
     232          if( (*aerprod) == 1 )
     233            {
     234              aer(corps,TEMP,NB,Y,&j,k_dep,faer,
     235                  &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
     236                  mmolaer,productaer,csurn,csurh);
     237
     238              for( i = 0; i <= 3; i++ )
     239                {
     240                  PRODAER[i][j] = productaer[i];
    230241                  MAER[i][j] = mmolaer[i];
    231                    CSN[i][j] = csurn[i];
    232                    CSH[i][j] = csurh[i];
    233              }
    234 /* DEBUG
    235 printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
    236 if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
    237       printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
    238               fp[utilaer[2]],dyc2h2*NB[j]);
    239 if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
    240       printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
    241               fp[utilaer[5]],dyhcn*NB[j]);
    242 if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
    243       printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
    244               fp[utilaer[6]],dyhc3n*NB[j]);
    245 if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
    246       printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
    247               fp[utilaer[13]],dynccn*NB[j]);
    248 if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
    249       printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
    250               fp[utilaer[14]],dych3cn*NB[j]);
    251 if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
    252       printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
    253               fp[utilaer[15]],dyc2h3cn*NB[j]);
    254 */
    255 
    256              fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
    257              fp[utilaer[5]] -= (    dyhcn * NB[j] );
    258              fp[utilaer[6]] -= (   dyhc3n * NB[j] );
    259              fp[utilaer[13]]-= (   dynccn * NB[j] );
    260              fp[utilaer[14]]-= (  dych3cn * NB[j] );
    261              fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
    262              if( Y[utilaer[2]][j]  != 0.0 )
    263        jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
    264              if( Y[utilaer[5]][j]  != 0.0 )
    265        jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
    266              if( Y[utilaer[6]][j]  != 0.0 )
    267        jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
    268              if( Y[utilaer[13]][j] != 0.0 )
    269      jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
    270              if( Y[utilaer[14]][j] != 0.0 )
    271      jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
    272              if( Y[utilaer[15]][j] != 0.0 )
    273      jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
    274          }
     242                  CSN[i][j] = csurn[i];
     243                  CSH[i][j] = csurh[i];
     244                }
     245              /* DEBUG
     246                 printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
     247                 if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
     248                printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
     249                fp[utilaer[2]],dyc2h2*NB[j]);
     250                 if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
     251                printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
     252                fp[utilaer[5]],dyhcn*NB[j]);
     253                 if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
     254                printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
     255                fp[utilaer[6]],dyhc3n*NB[j]);
     256                 if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
     257                printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
     258                fp[utilaer[13]],dynccn*NB[j]);
     259                 if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
     260                printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
     261                fp[utilaer[14]],dych3cn*NB[j]);
     262                 if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
     263                printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
     264                fp[utilaer[15]],dyc2h3cn*NB[j]);
     265              */
     266
     267              fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
     268              fp[utilaer[5]] -= (    dyhcn * NB[j] );
     269              fp[utilaer[6]] -= (   dyhc3n * NB[j] );
     270              fp[utilaer[13]]-= (   dynccn * NB[j] );
     271              fp[utilaer[14]]-= (  dych3cn * NB[j] );
     272              fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
     273              if( Y[utilaer[2]][j]  != 0.0 )
     274                jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
     275              if( Y[utilaer[5]][j]  != 0.0 )
     276                jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
     277              if( Y[utilaer[6]][j]  != 0.0 )
     278                jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
     279              if( Y[utilaer[13]][j] != 0.0 )
     280                jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
     281              if( Y[utilaer[14]][j] != 0.0 )
     282                jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
     283              if( Y[utilaer[15]][j] != 0.0 )
     284                jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
     285            }
    275286     
    276287       
    277 /* H -> H2 on haze particles */
    278 /* ------------------------- */
    279          if( (*htoh2) == 1 )
    280          {
     288          /* H -> H2 on haze particles */
     289          /* ------------------------- */
     290          if( (*htoh2) == 1 )
     291            {
    281292              heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
    282 /* dyh <= 0 / 1.0 en adsor., 1 en reac. */
    283 
    284 /* DEBUG
    285 printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
    286 if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
    287 printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
    288 if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
    289 printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
    290 */
     293              /* dyh <= 0 / 1.0 en adsor., 1 en reac. */
     294
     295              /* DEBUG
     296                 printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
     297                 if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
     298                 printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
     299                 if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
     300                 printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
     301              */
    291302
    292303              fp[utilaer[0]] += ( dyh  * NB[j] );
    293    /* pourquoi pas *2 ?? cf gptit dans 2da... */
     304              /* pourquoi pas *2 ?? cf gptit dans 2da... */
    294305
    295306              fp[utilaer[1]] += ( dyh2 * NB[j] );
    296307              if( Y[utilaer[0]][j] != 0.0 )
    297        jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
    298    /* pourquoi pas *2 ?? cf gptit dans 2da... */
    299          }
    300 
    301 
    302 /* Backup jacobian level j. */
    303 /* ------------------------ */
    304          for( i = 0; i <= ST-1; i++ )
     308                jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
     309              /* pourquoi pas *2 ?? cf gptit dans 2da... */
     310            }
     311
     312
     313          /* Backup jacobian level j. */
     314          /* ------------------------ */
     315          for( i = 0; i <= ST-1; i++ )
    305316            for( k = 0; k <= ST-1; k++ )
    306                a[j][i][k] = jac[i][k];     
    307 
    308 
    309 /* Diffusion verticale et flux exterieurs */
    310 /* -------------------------------------- */
    311 
    312 /*
    313 pour dy/dr, dr doit etre en cm...
    314 pareil pour dphi/dr
    315 */
    316          for( i = 0; i <= ST-1; i++ )
    317          {
    318 
    319 /* First level. */
    320             if( j == NLD )
    321             {
    322                v = dv = 0.0e0;
    323                dra = RA[j+1]-RA[j];
    324 
    325                cp  = (NB[j+1]+NB[j])/2.;  /* Mean total concentration. */
    326                dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
    327                      pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
    328                dp  = (MD[i][j]+MD[i][j+1])/2.;     /* Mean molecular diffusion. */
    329                dyp = (Y[i][j+1]-Y[i][j])/(RA[j+2]-RA[j])*2.e-5; /* Delta y level +1. */
    330                kp  = (KEDD[j+1]+KEDD[j])/2.;       /* Mean eddy diffusion. */
    331              /* div phi. */
    332                f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
    333                               + kp * dyp )
    334                         * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.))
    335                        + fp[i] - Y[i][j]*fl[i] + v;
    336              /* dphi / dy this level. */
    337                a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
    338                                     - 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
    339                         * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) + dv );
    340              /* dphi / dy level +1. */
    341                c[j][i]     = -THETA * delta
    342                              * cp * ( dp * 0.5e0 * dip
    343                                     + 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
    344                         * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.));
    345             }
    346 /* Last level. */
    347             else if( j == NLEV-1 )
    348             {
    349                v = dv = 0.0e0;
    350                dra = RA[NLEV-1]-RA[NLEV-2];
    351 
    352    /* Jeans escape */
    353                if( strcmp(corps[i], "H") == 0 ) 
    354                {
    355                  dv = top_H  * NB[NLEV-1]
     317              a[j][i][k] = jac[i][k];     
     318
     319
     320          /* Diffusion verticale et flux exterieurs */
     321          /* -------------------------------------- */
     322
     323          /*
     324            pour dy/dr, dr doit etre en cm...
     325            pareil pour dphi/dr
     326          */
     327          for( i = 0; i <= ST-1; i++ )
     328            {
     329
     330              /* First level. */
     331              if( j == NLD )
     332                {
     333                  v = dv = 0.0e0;
     334                  dra = RA[j+1]-RA[j];
     335
     336                  cp  = (NB[j+1]+NB[j])/2.;  /* Mean total concentration. */
     337                  dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
     338                    pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
     339                  dp  = (MD[i][j]+MD[i][j+1])/2.;     /* Mean molecular diffusion. */
     340                  dyp = (Y[i][j+1]-Y[i][j])/(RA[j+2]-RA[j])*2.e-5; /* Delta y level +1. */
     341                  kp  = (KEDD[j+1]+KEDD[j])/2.;       /* Mean eddy diffusion. */
     342                  /* div phi. */
     343                  f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
     344                                   + kp * dyp )
     345                    * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.))
     346                    + fp[i] - Y[i][j]*fl[i] + v;
     347                  /* dphi / dy this level. */
     348                  a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
     349                                        - 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
     350                                  * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.)) + dv );
     351                  /* dphi / dy level +1. */
     352                  c[j][i]     = -THETA * delta
     353                    * cp * ( dp * 0.5e0 * dip
     354                             + 2.e-5/(RA[j+2]-RA[j]) * (dp + kp) )
     355                    * (4.e-5/dra/pow((1.+RA[j]/RA[j+1]),2.));
     356                }
     357              /* Last level. */
     358              else if( j == NLEV-1 )
     359                {
     360                  v = dv = 0.0e0;
     361                  dra = RA[NLEV-1]-RA[NLEV-2];
     362
     363                  /* Jeans escape */
     364                  if( i == ih ) 
     365                    {
     366                      dv = top_H  * NB[NLEV-1]
    356367                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
    357                  v  = dv * Y[i][NLEV-1];
    358                }
    359                if( strcmp(corps[i], "H2") == 0 )
    360                {
    361                  dv = top_H2 * NB[NLEV-1]
     368                      v  = dv * Y[i][NLEV-1];
     369                    }
     370                  else if( i == ih2 )
     371                    {
     372                      dv = top_H2 * NB[NLEV-1]
    362373                        * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
    363                  v  = dv * Y[i][NLEV-1];
    364                }
    365    /* Input flux for N(4S) */
    366                if( strcmp(corps[i], "N4S") == 0 )
    367                  v  = top_N4S
    368                         * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
    369 
    370                cm  = (NB[NLEV-1]+NB[NLEV-2])/2.;  /* Mean total concentration. */
    371                dim = r * (MASS[i]-(mu[NLEV-1]+mu[NLEV-2])/2.)
    372                        / (TEMP[NLEV-1]+TEMP[NLEV-2])
    373                        / pow( RA[NLEV-1],   2.0e0 );  /* Delta i,j level -1. */
    374                dm  = (MD[i][NLEV-1]+MD[i][NLEV-2])/2.;    /* Mean molecular diffusion. */
    375                dym = (Y[i][NLEV-1]-Y[i][NLEV-2])/dra*1.e-5; /* Delta y level -1. */
    376                km  = (KEDD[NLEV-1]+KEDD[NLEV-2])/2.;      /* Mean eddy diffusion. */
    377              /* div phi. */
    378                f[i][NLEV-1] = fp[i] - Y[i][NLEV-1]*fl[i] - v
    379                        - cm * ( dm * ( (Y[i][NLEV-1]+Y[i][NLEV-2])/2. * dim + dym )
    380                               + km * dym )
    381                         * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
    382              /* dphi / dy this level */
    383                a[NLEV-1][i][i] -= ( cm * ( dm * 0.5e0 * dim
    384                                     + 1.e-5/dra * (dm + km ) )
    385                         * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)) + dv );
    386              /* dphi / dy level -1. */
    387                b[NLEV-1][i][2]  =  THETA * delta
    388                                   * cm * ( dm * 0.5e0 * dim
    389                                     - 1.e-5/dra * (dm + km ) )
    390                         * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
    391             }
    392             else
    393             {
    394                v = dv = 0.0e0;
    395                dram=(RA[j+1]-RA[j-1])/2.;
    396                if (j<NLEV-2)
    397                  drap=(RA[j+1]-RA[j-1])/2.;
    398                else
    399                  drap=dram;
    400 
    401                cm  = (NB[j]+NB[j-1])/2.;       /* Mean concentration level -1. */
    402                cp  = (NB[j]+NB[j+1])/2.;       /* Mean concentration level +1. */
    403                dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
    404                      pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
    405                dim = r * (MASS[i]-(mu[j]+mu[j-1])/2.) / (TEMP[j]+TEMP[j-1]) /
    406                      pow( RA[j],   2.0e0 );    /* Delta i,j level -1. */
    407                dm  = (MD[i][j-1]+MD[i][j])/2.;    /* Mean molecular diffusion level -1. */
    408                dp  = (MD[i][j+1]+MD[i][j])/2.;    /* Mean molecular diffusion level +1. */
    409                dym = (Y[i][j]-Y[i][j-1])/dram*1.e-5; /* Delta y level -1. */
    410                dyp = (Y[i][j+1]-Y[i][j])/drap*1.e-5; /* Delta y level +1. */
    411                km  = (KEDD[j]+KEDD[j-1])/2.;      /* Mean eddy diffusion level -1. */
    412                kp  = (KEDD[j]+KEDD[j+1])/2.;      /* Mean eddy diffusion level +1. */
    413              /* div phi. */
    414                f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
    415                               + kp * dyp )
    416                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
    417                        - cm * ( dm * ( (Y[i][j]+Y[i][j-1])/2. * dim + dym )
    418                               + km * dym )
    419                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.))
    420                        + fp[i] - fl[i] * Y[i][j] + v;
    421              /* dphi / dy this level */
    422                a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
    423                                     - 1.e-5/drap * (dp + kp) )
    424                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
    425                              - cm * ( dm * 0.5e0 * dim
    426                                     + 1.e-5/dram * (dm + km ) )
    427                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) );
    428              /* dphi / dy level -1. */
    429                b[j][i][2]  =  THETA * delta
    430                              * cm * ( dm * 0.5e0 * dim
    431                                     - 1.e-5/dram * (dm + km ) )
    432                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.));
    433              /* dphi / dy level +1. */
    434                c[j][i]     = -THETA * delta
    435                              * cp * ( dp * 0.5e0 * dip
    436                                     + 1.e-5/drap * (dp + kp) )
    437                         * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.));
    438             }
    439          }
    440 
    441 
    442 
    443 /* finition pour inversion */
    444 /* ----------------------- */
    445 
    446          for( i = 0; i <= ST-1; i++ )
    447          {
    448             for( k = 0; k <= ST-1; k++ )
    449             {
    450                a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
    451                if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
    452             }
    453             f[i][j] *= delta;
    454          }
    455 
    456    }
    457 
    458 
    459 /* -------------------------------- */
    460 /* Inversion of matrix cf method LU */
    461 /* -------------------------------- */
    462 
    463    for( j = NLD+1; j <= NLEV-1; j++ )
    464    {
    465          solve( a, j-1, 0, ST-1 );
    466          for( i = 0; i <= ST-1; i++ )
    467          {
    468             s = 0.0e0;
    469             for( k = 0; k <= ST-1; k++ )
    470             {
    471                a[j][i][k] -= ( b[j][i][2] * c[j-1][k] * a[j-1][i][k] );
    472                s          += ( b[j][i][2] * f[k][j-1] * a[j-1][i][k] );
    473             }
    474             f[i][j] -= s;
    475          }
    476    }
    477    solve( a, NLEV-1, 0, ST-1 );
    478    for( j = NLEV-1; j >= NLD; j-- )     
    479    {
    480          if( j != NLEV-1 )
     374                      v  = dv * Y[i][NLEV-1];
     375                    }
     376                  /* Input flux for N(4S) */
     377                  else if( i == in4s )
     378                    v  = top_N4S
     379                      * (4.e-5/dra/pow((2.-dra/(RA[NLEV-1]+dra)),2.));
     380
     381                  cm  = (NB[NLEV-1]+NB[NLEV-2])/2.;  /* Mean total concentration. */
     382                  dim = r * (MASS[i]-(mu[NLEV-1]+mu[NLEV-2])/2.)
     383                    / (TEMP[NLEV-1]+TEMP[NLEV-2])
     384                    / pow( RA[NLEV-1],   2.0e0 );  /* Delta i,j level -1. */
     385                  dm  = (MD[i][NLEV-1]+MD[i][NLEV-2])/2.;    /* Mean molecular diffusion. */
     386                  dym = (Y[i][NLEV-1]-Y[i][NLEV-2])/dra*1.e-5; /* Delta y level -1. */
     387                  km  = (KEDD[NLEV-1]+KEDD[NLEV-2])/2.;      /* Mean eddy diffusion. */
     388                  /* div phi. */
     389                  f[i][NLEV-1] = fp[i] - Y[i][NLEV-1]*fl[i] - v
     390                    - cm * ( dm * ( (Y[i][NLEV-1]+Y[i][NLEV-2])/2. * dim + dym )
     391                             + km * dym )
     392                    * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
     393                  /* dphi / dy this level */
     394                  a[NLEV-1][i][i] -= ( cm * ( dm * 0.5e0 * dim
     395                                              + 1.e-5/dra * (dm + km ) )
     396                                       * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.)) + dv );
     397                  /* dphi / dy level -1. */
     398                  b[NLEV-1][i][2]  =  THETA * delta
     399                    * cm * ( dm * 0.5e0 * dim
     400                             - 1.e-5/dra * (dm + km ) )
     401                    * (4.e-5/dra/pow((2.+dra/RA[NLEV-1]),2.));
     402                }
     403              else
     404                {
     405                  v = dv = 0.0e0;
     406                  dram=(RA[j+1]-RA[j-1])/2.;
     407                  if (j<NLEV-2)
     408                    drap=(RA[j+1]-RA[j-1])/2.;
     409                  else
     410                    drap=dram;
     411
     412                  cm  = (NB[j]+NB[j-1])/2.;       /* Mean concentration level -1. */
     413                  cp  = (NB[j]+NB[j+1])/2.;       /* Mean concentration level +1. */
     414                  dip = r * (MASS[i]-(mu[j+1]+mu[j])/2.) / (TEMP[j+1]+TEMP[j]) /
     415                    pow( RA[j+1], 2.0e0 );    /* Delta i,j level +1. */
     416                  dim = r * (MASS[i]-(mu[j]+mu[j-1])/2.) / (TEMP[j]+TEMP[j-1]) /
     417                    pow( RA[j],   2.0e0 );    /* Delta i,j level -1. */
     418                  dm  = (MD[i][j-1]+MD[i][j])/2.;    /* Mean molecular diffusion level -1. */
     419                  dp  = (MD[i][j+1]+MD[i][j])/2.;    /* Mean molecular diffusion level +1. */
     420                  dym = (Y[i][j]-Y[i][j-1])/dram*1.e-5; /* Delta y level -1. */
     421                  dyp = (Y[i][j+1]-Y[i][j])/drap*1.e-5; /* Delta y level +1. */
     422                  km  = (KEDD[j]+KEDD[j-1])/2.;      /* Mean eddy diffusion level -1. */
     423                  kp  = (KEDD[j]+KEDD[j+1])/2.;      /* Mean eddy diffusion level +1. */
     424                  /* div phi. */
     425                  f[i][j] = cp * ( dp * ( (Y[i][j+1]+Y[i][j])/2. * dip + dyp )
     426                                   + kp * dyp )
     427                    * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
     428                    - cm * ( dm * ( (Y[i][j]+Y[i][j-1])/2. * dim + dym )
     429                             + km * dym )
     430                    * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.))
     431                    + fp[i] - fl[i] * Y[i][j] + v;
     432                  /* dphi / dy this level */
     433                  a[j][i][i] += ( cp * ( dp * 0.5e0 * dip
     434                                        - 1.e-5/drap * (dp + kp) )
     435                                  * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.))
     436                                  - cm * ( dm * 0.5e0 * dim
     437                                           + 1.e-5/dram * (dm + km ) )
     438                                  * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.)) );
     439                  /* dphi / dy level -1. */
     440                  b[j][i][2]  =  THETA * delta
     441                    * cm * ( dm * 0.5e0 * dim
     442                             - 1.e-5/dram * (dm + km ) )
     443                    * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j+1]/RA[j]),2.));
     444                  /* dphi / dy level +1. */
     445                  c[j][i]     = -THETA * delta
     446                    * cp * ( dp * 0.5e0 * dip
     447                             + 1.e-5/drap * (dp + kp) )
     448                    * (4.e-5/(RA[j+1]-RA[j])/pow((1.+RA[j]/RA[j+1]),2.));
     449                }
     450            }
     451
     452
     453
     454          /* finition pour inversion */
     455          /* ----------------------- */
     456
     457          for( i = 0; i <= ST-1; i++ )
     458            {
     459              for( k = 0; k <= ST-1; k++ )
     460                {
     461                  a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
     462                  if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
     463                }
     464              f[i][j] *= delta;
     465            }
     466
     467        }
     468
     469
     470      /* -------------------------------- */
     471      /* Inversion of matrix cf method LU */
     472      /* -------------------------------- */
     473
     474      for( j = NLD+1; j <= NLEV-1; j++ )
     475        {
     476          solve( a, j-1, 0, ST-1 );
     477          for( i = 0; i <= ST-1; i++ )
     478            {
     479              s = 0.0e0;
     480              for( k = 0; k <= ST-1; k++ )
     481                {
     482                  a[j][i][k] -= ( b[j][i][2] * c[j-1][k] * a[j-1][i][k] );
     483                  s          += ( b[j][i][2] * f[k][j-1] * a[j-1][i][k] );
     484                }
     485              f[i][j] -= s;
     486            }
     487        }
     488      solve( a, NLEV-1, 0, ST-1 );
     489      for( j = NLEV-1; j >= NLD; j-- )     
     490        {
     491          if( j != NLEV-1 )
    481492            for( i = 0; i <= ST-1; i++ ) f[i][j] -= ( c[j][i] * b[j+1][i][1] );
    482          for( i = 0; i <= ST-1; i++ )
    483          {
    484             s = 0.0e0;
    485             for( k = 0; k <= ST-1; k++ ) s += ( a[j][i][k] * f[k][j] );
    486             b[j][i][1]  = s;
    487             Y[i][j]    += s;
    488             if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
    489          }
    490    }
    491 
    492 /* ------------------ */
    493 /* Tests et evolution */
    494 /* ------------------ */
    495 
    496 /* Calcul deviation */
    497 /* ---------------- */
    498 
    499    for( j = NLD; j <= NLEV-1; j++ )
    500       for( i = 0; i <= ST-1; i++ )
    501          if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
    502          {
    503                conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
    504                if( conv > ts )
    505                {
    506 /*
    507                   if( conv >= 0.1 )
    508                   {
    509                      out = fopen( outlog, "a" );
    510                      fprintf( out, "Latitude %g;", (*LAT));
    511                      fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
    512                      fclose( out );
    513                   }
    514 */
     493          for( i = 0; i <= ST-1; i++ )
     494            {
     495              s = 0.0e0;
     496              for( k = 0; k <= ST-1; k++ ) s += ( a[j][i][k] * f[k][j] );
     497              b[j][i][1]  = s;
     498              Y[i][j]    += s;
     499              if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
     500            }
     501        }
     502
     503      /* ------------------ */
     504      /* Tests et evolution */
     505      /* ------------------ */
     506
     507      /* Calcul deviation */
     508      /* ---------------- */
     509
     510      for( j = NLD; j <= NLEV-1; j++ )
     511        for( i = 0; i <= ST-1; i++ )
     512          if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
     513            {
     514              conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
     515              if( conv > ts )
     516                {
     517                  /*
     518                    if( conv >= 0.1 )
     519                    {
     520                    out = fopen( outlog, "a" );
     521                    fprintf( out, "Latitude %g;", (*LAT));
     522                    fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
     523                    fclose( out );
     524                    }
     525                  */
    515526                  ts = conv;
    516                }
    517          }
    518 
    519 /* test deviation */
    520 /* -------------- */
    521 
    522          if( ts < 0.1e0 )
    523          {
    524             for( i = 0; i <= ST-1; i++ )
    525                for( j = NLD; j <= NLEV-1; j++ )
    526                  if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
    527                  {
     527                }
     528            }
     529
     530      /* test deviation */
     531      /* -------------- */
     532
     533      if( ts < 0.1e0 )
     534        {
     535          for( i = 0; i <= ST-1; i++ )
     536            for( j = NLD; j <= NLEV-1; j++ )
     537              if( (Y[i][j] >= 0.5e0) && (i != in2) )
     538                {
    528539                  out = fopen( outlog, "a" );
    529540                  fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
    … …  
    532543                  fclose( out );
    533544                  exit(0);
    534 //                  Y[i][j] = 1.e-20;
    535                  }
    536             for( j = NLD; j <= NLEV-1; j++ )
    537                for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
    538             time += delta;
    539             if(   ts < 1.00e-5 )                      delta *= 1.0e2;
    540             if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
    541             if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
    542             if( ( ts > 1.00e-3 ) && ( ts < 5.0e-3 ) ) delta *= 3.0e0;
    543             if( ( ts > 5.00e-3 ) && ( ts < 0.01e0 ) ) delta *= 1.5e0;
    544             if( ( ts > 0.010e0 ) && ( ts < 0.03e0 ) ) delta *= 1.2e0;
    545             if( ( ts > 0.030e0 ) && ( ts < 0.05e0 ) ) delta *= 1.1e0;
    546 
    547 //            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
    548 //            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
     545                  //                  Y[i][j] = 1.e-20;
     546                }
     547          for( j = NLD; j <= NLEV-1; j++ )
     548            for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
     549          time += delta;
     550          if(   ts < 1.00e-5 )                      delta *= 1.0e2;
     551          if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
     552          if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
     553          if( ( ts > 1.00e-3 ) && ( ts < 5.0e-3 ) ) delta *= 3.0e0;
     554          if( ( ts > 5.00e-3 ) && ( ts < 0.01e0 ) ) delta *= 1.5e0;
     555          if( ( ts > 0.010e0 ) && ( ts < 0.03e0 ) ) delta *= 1.2e0;
     556          if( ( ts > 0.030e0 ) && ( ts < 0.05e0 ) ) delta *= 1.1e0;
     557
     558          //            if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
     559          //            if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
    549560         
    550             delta = min( deltamax, delta );
    551          }
    552          else
    553          {
    554             for( j = NLD; j <= NLEV-1; j++ )
    555                for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
    556 
    557             if(   ts > 0.8 )                    delta *= 1.e-6;
    558             if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
    559             if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
    560             if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
    561             if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
    562             if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
    563          }
    564          ts = 0.0e0;
    565 
    566          out = fopen( outlog, "a" );
    567          fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
    568          fclose( out );
    569 
    570 }
    571 /* **************** */       
    572 /* end of time loop */
    573 /* **************** */       
    574 
    575 /*
    576 ==========================================================================
    577 
    578 SECONDE ETAPE:
    579 ===============
     561          delta = min( deltamax, delta );
     562        }
     563      else
     564        {
     565          for( j = NLD; j <= NLEV-1; j++ )
     566            for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
     567
     568          if(   ts > 0.8 )                    delta *= 1.e-6;
     569          if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
     570          if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
     571          if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
     572          if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
     573          if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
     574        }
     575      ts = 0.0e0;
     576
     577      out = fopen( outlog, "a" );
     578      fprintf(out, "delta:%e; time:%e; fin:%e\n",delta,time,(*FIN));
     579      fclose( out );
     580
     581    }
     582  /* **************** */       
     583  /* end of time loop */
     584  /* **************** */       
     585
     586  /*
     587    ==========================================================================
     588
     589    SECONDE ETAPE:
     590    ===============
    580591    INVERSION LOCALE PAR BLOC ENTRE NLD ET LA SURFACE
    581 ===============
    582 */
    583    if( NLD != 0 )
    584    for( j = NLD-1; j >= 0; j-- )
    585    {
    586       time     = ts = 0.0e0;
    587       delta    = 1.e-3;
    588 
    589 /* ++++++++++++ */
    590 /*  time loop.  */
    591 /* ++++++++++++ */
    592 
    593       while( time < (*FIN) )     
     592    ===============
     593  */
     594  if( NLD != 0 )
     595    for( j = NLD-1; j >= 0; j-- )
    594596      {
    595 
    596 /* init of step */
    597 /* ------------ */
    598          for( i = 0; i <= ST-1; i++ )
    599          {
    600             fp[i] = fl[i] = 0.0e0;
    601             for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
    602          }
    603 
    604 /* Chimie */
    605 /* ------ */
    606 
    607 /* productions et pertes chimiques */
    608          for( i = 0; i <= ST-1; i++ )   
    609          {
    610             Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
    611 
    612             for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
    613             {
    614                ireac = prod[i][l];                  /* Number of the reaction involves. */
    615                ncom1 = reactif[ireac][0];           /* First compound which reacts. */
    616                if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
    617                {
    618                   jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
    619                   fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
    620                }
    621                else                                 /* General case. */
    622                {
    623                   ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
    624                   jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
    625                   jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
    626                   fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
    627                }
    628             }
     597        time     = ts = 0.0e0;
     598        delta    = 1.e-3;
     599
     600        /* ++++++++++++ */
     601        /*  time loop.  */
     602        /* ++++++++++++ */
     603
     604        while( time < (*FIN) )     
     605          {
     606
     607            /* init of step */
     608            /* ------------ */
     609            for( i = 0; i <= ST-1; i++ )
     610              {
     611                fp[i] = fl[i] = 0.0e0;
     612                for( l = 0; l <= ST-1; l++ ) jac[i][l] = 0.0e0;
     613              }
     614
     615            /* Chimie */
     616            /* ------ */
     617
     618            /* productions et pertes chimiques */
     619            for( i = 0; i <= ST-1; i++ )   
     620              {
     621                Y[i][j] = max(Y[i][j],1.e-30);                /* minimum */
     622
     623                for( l = 0; l <= nom_prod[i]-1; l++ )    /* Production term */
     624                  {
     625                    ireac = prod[i][l];                  /* Number of the reaction involves. */
     626                    ncom1 = reactif[ireac][0];           /* First compound which reacts. */
     627                    if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     628                      {
     629                        jac[i][ncom1] += ( KRATE[ireac][j] * NB[j] );
     630                        fp[i]         += ( KRATE[ireac][j] * NB[j] * Y[ncom1][j] );
     631                      }
     632                    else                                 /* General case. */
     633                      {
     634                        ncom2          = reactif[ireac][1];                       /* Second compound which reacts. */
     635                        jac[i][ncom1] += ( KRATE[ireac][j] * Y[ncom2][j] );       /* Jacobian compound #1. */
     636                        jac[i][ncom2] += ( KRATE[ireac][j] * Y[ncom1][j] );       /* Jacobian compound #2. */
     637                        fp[i] += ( KRATE[ireac][j] * Y[ncom1][j] * Y[ncom2][j] ); /* Production term. */
     638                      }
     639                  }
    629640           
    630             for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
    631             {
    632                ireac = perte[i][l][0];              /* Reaction number. */
    633                ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
    634                if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
    635                {
    636                   jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
    637                   fl[i]     += ( KRATE[ireac][j] * NB[j] );
    638                }
    639                else                                 /* General case. */
    640                {
    641                   jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
    642                   jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
    643                   fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
    644                }
    645             }
    646          }
    647 
    648 
    649 /* Aerosols */
    650 /* -------- */
    651          if( (*aerprod) == 1 )
    652          {
    653              aer(corps,TEMP,NB,Y,&j,k_dep,faer,
    654               &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
    655               mmolaer,productaer,csurn,csurh);
    656 
    657              for( i = 0; i <= 3; i++ )
    658              {
    659                PRODAER[i][j] = productaer[i];
    660                   MAER[i][j] = mmolaer[i];
    661                    CSN[i][j] = csurn[i];
    662                    CSH[i][j] = csurh[i];
    663              }
    664 /* DEBUG
    665 printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
    666 if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
    667       printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
    668               fp[utilaer[2]],dyc2h2*NB[j]);
    669 if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
    670       printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
    671               fp[utilaer[5]],dyhcn*NB[j]);
    672 if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
    673       printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
    674               fp[utilaer[6]],dyhc3n*NB[j]);
    675 if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
    676       printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
    677               fp[utilaer[13]],dynccn*NB[j]);
    678 if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
    679       printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
    680               fp[utilaer[14]],dych3cn*NB[j]);
    681 if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
    682       printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
    683               fp[utilaer[15]],dyc2h3cn*NB[j]);
    684 */
    685 
    686              fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
    687              fp[utilaer[5]] -= (    dyhcn * NB[j] );
    688              fp[utilaer[6]] -= (   dyhc3n * NB[j] );
    689              fp[utilaer[13]]-= (   dynccn * NB[j] );
    690              fp[utilaer[14]]-= (  dych3cn * NB[j] );
    691              fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
    692              if( Y[utilaer[2]][j]  != 0.0 )
    693        jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
    694              if( Y[utilaer[5]][j]  != 0.0 )
    695        jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
    696              if( Y[utilaer[6]][j]  != 0.0 )
    697        jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
    698              if( Y[utilaer[13]][j] != 0.0 )
    699      jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
    700              if( Y[utilaer[14]][j] != 0.0 )
    701      jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
    702              if( Y[utilaer[15]][j] != 0.0 )
    703      jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
    704          }
     641                for( l = 0; l <= nom_perte[i]-1; l++ )   /* Loss term. */
     642                  {
     643                    ireac = perte[i][l][0];              /* Reaction number. */
     644                    ncom2 = perte[i][l][1];              /* Compound #2 reacts. */
     645                    if( reactif[ireac][1] == NC )        /* Photodissociation or relaxation */
     646                      {
     647                        jac[i][i] -= ( KRATE[ireac][j] * NB[j] );
     648                        fl[i]     += ( KRATE[ireac][j] * NB[j] );
     649                      }
     650                    else                                 /* General case. */
     651                      {
     652                        jac[i][ncom2] -= ( KRATE[ireac][j] * Y[i][j] );       /* Jacobian compound #1. */
     653                        jac[i][i]     -= ( KRATE[ireac][j] * Y[ncom2][j] );   /* Jacobien compound #2. */
     654                        fl[i]         += ( KRATE[ireac][j] * Y[ncom2][j] );   /* Loss term. */
     655                      }
     656                  }
     657              }
     658
     659
     660            /* Aerosols */
     661            /* -------- */
     662            if( (*aerprod) == 1 )
     663              {
     664                aer(corps,TEMP,NB,Y,&j,k_dep,faer,
     665                    &dyc2h2,&dyhc3n,&dyhcn,&dynccn,&dych3cn,&dyc2h3cn,utilaer,
     666                    mmolaer,productaer,csurn,csurh);
     667
     668                for( i = 0; i <= 3; i++ )
     669                  {
     670                    PRODAER[i][j] = productaer[i];
     671                    MAER[i][j] = mmolaer[i];
     672                    CSN[i][j] = csurn[i];
     673                    CSH[i][j] = csurh[i];
     674                  }
     675                /* DEBUG
     676                   printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
     677                   if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
     678                   printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
     679                   fp[utilaer[2]],dyc2h2*NB[j]);
     680                   if(fabs(dyhcn*NB[j])>fabs(fp[utilaer[5]]/10.))
     681                   printf("fp(%s) =%e; dyhcn  =%e\n",corps[utilaer[5]],
     682                   fp[utilaer[5]],dyhcn*NB[j]);
     683                   if(fabs(dyhc3n*NB[j])>fabs(fp[utilaer[6]]/10.))
     684                   printf("fp(%s) =%e; dyhc3n =%e\n",corps[utilaer[6]],
     685                   fp[utilaer[6]],dyhc3n*NB[j]);
     686                   if(fabs(dynccn*NB[j])>fabs(fp[utilaer[13]]/10.))
     687                   printf("fp(%s) =%e; dynccn =%e\n",corps[utilaer[13]],
     688                   fp[utilaer[13]],dynccn*NB[j]);
     689                   if(fabs(dych3cn*NB[j])>fabs(fp[utilaer[14]]/10.))
     690                   printf("fp(%s) =%e; dych3cn=%e\n",corps[utilaer[14]],
     691                   fp[utilaer[14]],dych3cn*NB[j]);
     692                   if(fabs(dyc2h3cn*NB[j])>fabs(fp[utilaer[15]]/10.))
     693                   printf("fp(%s) =%e; dyc2h3cn=%e\n",corps[utilaer[15]],
     694                   fp[utilaer[15]],dyc2h3cn*NB[j]);
     695                */
     696
     697                fp[utilaer[2]] -= (   dyc2h2 * NB[j] );
     698                fp[utilaer[5]] -= (    dyhcn * NB[j] );
     699                fp[utilaer[6]] -= (   dyhc3n * NB[j] );
     700                fp[utilaer[13]]-= (   dynccn * NB[j] );
     701                fp[utilaer[14]]-= (  dych3cn * NB[j] );
     702                fp[utilaer[15]]-= ( dyc2h3cn * NB[j] );
     703                if( Y[utilaer[2]][j]  != 0.0 )
     704                  jac[utilaer[2]][utilaer[2]] -= (  dyc2h2 * NB[j] / Y[utilaer[2]][j] );
     705                if( Y[utilaer[5]][j]  != 0.0 )
     706                  jac[utilaer[5]][utilaer[5]] -= (   dyhcn * NB[j] / Y[utilaer[5]][j] );
     707                if( Y[utilaer[6]][j]  != 0.0 )
     708                  jac[utilaer[6]][utilaer[6]] -= (  dyhc3n * NB[j] / Y[utilaer[6]][j] );
     709                if( Y[utilaer[13]][j] != 0.0 )
     710                  jac[utilaer[13]][utilaer[13]] -= (  dynccn * NB[j] / Y[utilaer[13]][j] );
     711                if( Y[utilaer[14]][j] != 0.0 )
     712                  jac[utilaer[14]][utilaer[14]] -= ( dych3cn * NB[j] / Y[utilaer[14]][j] );
     713                if( Y[utilaer[15]][j] != 0.0 )
     714                  jac[utilaer[15]][utilaer[15]] -= (dyc2h3cn * NB[j] / Y[utilaer[15]][j] );
     715              }
    705716     
    706717       
    707 /* H -> H2 on haze particles */
    708 /* ------------------------- */
    709          if( (*htoh2) == 1 )
    710          {
    711               heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
    712 /* dyh <= 0 / 1.0 en adsor., 1 en reac. */
    713 
    714 /* DEBUG
    715 printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
    716 if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
    717 printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
    718 if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
    719 printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
    720 */
    721 
    722               fp[utilaer[0]] += ( dyh  * NB[j] );
    723    /* pourquoi pas *2 ?? cf gptit dans 2da... */
    724 
    725               fp[utilaer[1]] += ( dyh2 * NB[j] );
    726               if( Y[utilaer[0]][j] != 0.0 )
    727        jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
    728    /* pourquoi pas *2 ?? cf gptit dans 2da... */
    729          }
    730 
    731 
    732 /* Backup jacobian level j. */
    733 /* ------------------------ */
    734          for( i = 0; i <= ST-1; i++ )
    735          {
    736             for( k = 0; k <= ST-1; k++ )
    737                a[j][i][k] = jac[i][k];     
    738             f[i][j] = fp[i] - fl[i] * Y[i][j];
    739          }
    740 
    741 
    742 /* finition pour inversion */
    743 /* ----------------------- */
    744 
    745          for( i = 0; i <= ST-1; i++ )
    746          {
    747             for( k = 0; k <= ST-1; k++ )
    748             {
    749                a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
    750                if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
    751             }
    752             f[i][j] *= delta;
    753          }
    754 
    755 
    756 /* Inversion of matrix cf method LU */
    757 /* -------------------------------- */
    758 
    759 /* inversion by blocs: */
    760 /* Hydrocarbons */
    761 
    762          solve_b( a, f, j, 0, NHC-1 );             
    763          for( i = 0; i <= NHC-1; i++ )
    764          {
    765             Y[i][j] += f[i][j];
    766             if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
    767          }
    768 
    769 /* Nitriles */
    770 
    771          solve_b( a, f, j, NHC, ST-1 );             
    772          for( i = NHC+1; i <= ST-1; i++ )
    773          {
    774             Y[i][j] += f[i][j];
    775             if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
    776          }
    777 
    778 /* end inversion by blocs: */
    779 
    780 /* CH4 au sol */
    781 /* ---------- */
    782    for( i = 0; i <= ST-1; i++ )
    783      if( ( strcmp(corps[i], "CH4") == 0 ) && (j==0) && ( Y[i][0] < *botCH4 ) )
    784      {
    785           fluxCH4 += (*botCH4 - Y[i][0]);
    786            Y[i][0] = *botCH4;
    787      }
    788 
    789 /* ------------------ */
    790 /* Tests et evolution */
    791 /* ------------------ */
    792 
    793 /* Calcul deviation */
    794 /* ---------------- */
    795 
    796          for( i = 0; i <= ST-1; i++ )
    797          {
    798             test = 1.0e-15;
    799             if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
    800             {
    801                conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
    802 
    803                if( conv > ts )
    804                {
    805 /*
    806                   if( conv >= 0.1 )
    807                   {
    808                      out = fopen( outlog, "a" );
    809                      fprintf( out, "Latitude: %g; declin:%e;", (*LAT), (*DECLIN) );
    810                      fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
    811                      fclose( out );
    812                   }
    813 */
    814                   ts = conv;
    815                }
    816             }
    817          }
    818 
    819 /* test deviation */
    820 /* -------------- */
    821 
    822          if( ts < 0.1e0 )
    823          {
    824             for( i = 0; i <= ST-1; i++ )
    825                  if( (Y[i][j] >= 0.5e0) && (strcmp(corps[i],"N2") != 0) )
    826                  {
    827                   out = fopen( outlog, "a" );
    828                   fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
    829                            corps[i], ym1[i][j], Y[i][j], j );
    830                   for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i][j],Y[i][k] );
    831                   fclose( out );
    832              //     exit(0);
    833                   Y[i][j] = 1.e-20;
    834                  }
    835             for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
    836             time += delta;
    837             if(   ts < 1.00e-5 )                      delta *= 1.0e2;
    838             if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
    839             if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
    840             if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
    841             if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
     718            /* H -> H2 on haze particles */
     719            /* ------------------------- */
     720            if( (*htoh2) == 1 )
     721              {
     722                heterohtoh2(corps,TEMP,NB,Y,surfhaze,&j,&dyh,&dyh2,utilaer);
     723                /* dyh <= 0 / 1.0 en adsor., 1 en reac. */
     724
     725                /* DEBUG
     726                   printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
     727                   if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
     728                   printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
     729                   if(fabs(dyh2*NB[j])>fabs(fp[utilaer[1]]/10.))
     730                   printf("fp(%s) = %e; dyh2 = %e\n",corps[utilaer[1]],fp[utilaer[1]],dyh2*NB[j]);
     731                */
     732
     733                fp[utilaer[0]] += ( dyh  * NB[j] );
     734                /* pourquoi pas *2 ?? cf gptit dans 2da... */
     735
     736                fp[utilaer[1]] += ( dyh2 * NB[j] );
     737                if( Y[utilaer[0]][j] != 0.0 )
     738                  jac[utilaer[0]][utilaer[0]] += ( dyh  * NB[j] / Y[utilaer[0]][j] );
     739                /* pourquoi pas *2 ?? cf gptit dans 2da... */
     740              }
     741
     742
     743            /* Backup jacobian level j. */
     744            /* ------------------------ */
     745            for( i = 0; i <= ST-1; i++ )
     746              {
     747                for( k = 0; k <= ST-1; k++ )
     748                  a[j][i][k] = jac[i][k];     
     749                f[i][j] = fp[i] - fl[i] * Y[i][j];
     750              }
     751
     752
     753            /* finition pour inversion */
     754            /* ----------------------- */
     755
     756            for( i = 0; i <= ST-1; i++ )
     757              {
     758                for( k = 0; k <= ST-1; k++ )
     759                  {
     760                    a[j][i][k] *= ( -THETA * delta );  /* Correction time step. */
     761                    if( k == i ) a[j][k][k] += NB[j];  /* Correction diagonal. */
     762                  }
     763                f[i][j] *= delta;
     764              }
     765
     766
     767            /* Inversion of matrix cf method LU */
     768            /* -------------------------------- */
     769
     770            /* inversion by blocs: */
     771            /* Hydrocarbons */
     772
     773            solve_b( a, f, j, 0, NHC-1 );             
     774            for( i = 0; i <= NHC-1; i++ )
     775              {
     776                Y[i][j] += f[i][j];
     777                if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
     778              }
     779
     780            /* Nitriles */
     781
     782            solve_b( a, f, j, NHC, ST-1 );             
     783            for( i = NHC+1; i <= ST-1; i++ )
     784              {
     785                Y[i][j] += f[i][j];
     786                if( Y[i][j] <= 1.0e-30 ) Y[i][j] = 0.0e0;
     787              }
     788
     789            /* end inversion by blocs: */
     790
     791            /* CH4 au sol */
     792            /* ---------- */
     793            if( Y[ich4][0] < *botCH4 )
     794                {
     795                  dyCH4surf += (*botCH4 - Y[ich4][0]*10.0e0);
     796                  Y[ich4][0] = *botCH4;
     797                }
     798
     799            /* ------------------ */
     800            /* Tests et evolution */
     801            /* ------------------ */
     802
     803            /* Calcul deviation */
     804            /* ---------------- */
     805
     806            for( i = 0; i <= ST-1; i++ )
     807              {
     808                test = 1.0e-15;
     809                if( ( Y[i][j] > test ) && ( ym1[i][j] > test ) )
     810                  {
     811                    conv = fabs( Y[i][j] - ym1[i][j] ) / ym1[i][j];
     812
     813                    if( conv > ts )
     814                      {
     815                        /*
     816                          if( conv >= 0.1 )
     817                          {
     818                          out = fopen( outlog, "a" );
     819                          fprintf( out, "Latitude: %g; declin:%e;", (*LAT), (*DECLIN) );
     820                          fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
     821                          fclose( out );
     822                          }
     823                        */
     824                        ts = conv;
     825                      }
     826                  }
     827              }
     828
     829            /* test deviation */
     830            /* -------------- */
     831
     832            if( ts < 0.1e0 )
     833              {
     834                for( i = 0; i <= ST-1; i++ )
     835                  if( (Y[i][j] >= 0.5e0) && (i != in2) )
     836                    {
     837                      out = fopen( outlog, "a" );
     838                      fprintf( out, "WARNING %s mixing ratio is %e %e at %d\n",
     839                               corps[i], ym1[i][j], Y[i][j], j );
     840                      for( k = 0; k <= NLEV-1; k++ ) fprintf( out, "%d %e %e\n",k,ym1[i][j],Y[i][k] );
     841                      fclose( out );
     842                      //     exit(0);
     843                      Y[i][j] = 1.e-20;
     844                    }
     845                for( i = 0; i <= NC-1; i++ ) ym1[i][j] = max(Y[i][j],1.e-30);
     846                time += delta;
     847                if     (   ts < 1.00e-5 )                      delta *= 1.0e2;
     848                else if( ( ts > 1.00e-5 ) && ( ts < 1.0e-4 ) ) delta *= 1.0e1;
     849                else if( ( ts > 1.00e-4 ) && ( ts < 1.0e-3 ) ) delta *= 5.0e0;
     850                else if( ( ts > 0.001e0 ) && ( ts < 0.01e0 ) ) delta *= 3.0e0;
     851                else if( ( ts > 0.010e0 ) && ( ts < 0.05e0 ) ) delta *= 1.5e0;
    842852         
    843             delta = min( deltamax, delta );
    844          }
    845          else
    846          {
    847             for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
    848 
    849             if(   ts > 0.8 )                    delta *= 1.e-6;
    850             if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
    851             if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
    852             if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
    853             if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
    854             if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
    855          }
    856          ts = 0.0e0;
    857 /*
    858                      out = fopen( outlog, "a" );
    859                      fprintf(out, " alt:%e; delta:%e; time:%e; fin:%e\n",(RA[j]-R0),delta,time,(*FIN));
    860                      fclose( out );
    861 */
    862       }               
    863 
    864 /* +++++++++++++++++++ */
    865 /*  end of time loop.  */
    866 /* +++++++++++++++++++ */
    867 
    868       for( i = 0; i <= ST-1; i++ )
    869          if( ( strcmp(corps[i],"CH4") == 0 ) && ( j == 0 ) )
    870             fluxCH4 *= ( MASS[i]/(6.022e23*time) );
    871 
    872    }  /*  boucle j */
    873 
    874 
    875 /*
    876 ==========================================================================
    877 
    878 FINALISATION:
    879 ===============
    880 */
    881       for( i = 0; i <= ST-1; i++ )
    882          if( strcmp(corps[i],"CH4") == 0 )
    883             fluxCH4 *= ( MASS[i]/(6.022e23*time) );
    884 
    885 /* Niveau de N2 */
    886 /* ------------ */
     853                delta = min( deltamax, delta );
     854              }
     855            else
     856              {
     857                for( i = 0; i <= NC-1; i++ ) Y[i][j] = ym1[i][j];
     858
     859                if     (   ts > 0.8 )                    delta *= 1.e-6;
     860                else if( ( ts > 0.6 ) && ( ts <= 0.8 ) ) delta *= 1.e-4;
     861                else if( ( ts > 0.4 ) && ( ts <= 0.6 ) ) delta *= 1.e-2;
     862                else if( ( ts > 0.3 ) && ( ts <= 0.4 ) ) delta *= 0.1;
     863                else if( ( ts > 0.2 ) && ( ts <= 0.3 ) ) delta *= 0.2;
     864                else if( ( ts > 0.1 ) && ( ts <= 0.2 ) ) delta *= 0.3;
     865              }
     866            ts = 0.0e0;
     867            /*
     868              out = fopen( outlog, "a" );
     869              fprintf(out, " alt:%e; delta:%e; time:%e; fin:%e\n",(RA[j]-R0),delta,time,(*FIN));
     870              fclose( out );
     871            */
     872          }               
     873
     874        /* +++++++++++++++++++ */
     875        /*  end of time loop.  */
     876        /* +++++++++++++++++++ */
     877
     878      }  /*  boucle j */
     879
     880
     881  /*
     882    ==========================================================================
     883
     884    FINALISATION:
     885    ===============
     886
     887  /* Niveau de N2 */
     888  /* ------------ */
     889
     890  *fluxCH4=dyCH4surf;
    887891     
    888    for( j = 0; j <= NLEV-1; j++ )
    889    {
     892  for( j = 0; j <= NLEV-1; j++ )
     893    {
    890894      conv = 0.0e0;
    891895      for( i = 0; i <= ST-1; i++ )
    892          if( strcmp(corps[i],"N2") != 0 ) conv += Y[i][j];
    893       for( i = 0; i <= ST-1; i++ )
    894          if( strcmp(corps[i],"N2") == 0 ) Y[i][j] = 1. - conv;
    895    }
    896 
    897    if( (*aerprod) == 1 )
    898    {
    899     fdm2d( k_dep, 1, 5, 1 );
    900     fdm2d(  faer, 1, 5, 1 );
    901     fdm1d( productaer,  0 );
    902     fdm1d( mmolaer,  0 );
    903     fdm1d( csurn, 0 );
    904     fdm1d( csurh, 0 );
    905    }
    906 
    907    fdm1d(      fl, 0 );
    908    fdm1d(      fp, 0 );
    909    fdm1d(      mu, 0 );
    910    fdm2d(     ym1, 0,   NC-1, 0 );
    911    fdm2d(       f, 0,   NC-1, 0 );
    912    fdm2d(     jac, 0,   NC-1, 0 );
    913    fdm2d(       c, 0, NLEV-1, 0 );
    914    fdm3d(       a, 0, NLEV-1, 0,   NC-1, 0 );
    915    fdm3d(       b, 0, NLEV-1, 0,   NC-1, 1 );
     896        if( i != in2 ) conv += Y[i][j];
     897      Y[in2][j] = 1. - conv;
     898    }
     899
     900  if( (*aerprod) == 1 )
     901    {
     902      fdm2d( k_dep, 1, 5, 1 );
     903      fdm2d(  faer, 1, 5, 1 );
     904      fdm1d( productaer,  0 );
     905      fdm1d( mmolaer,  0 );
     906      fdm1d( csurn, 0 );
     907      fdm1d( csurh, 0 );
     908    }
     909
     910  fdm1d(      fl, 0 );
     911  fdm1d(      fp, 0 );
     912  fdm1d(      mu, 0 );
     913  fdm2d(     ym1, 0,   NC-1, 0 );
     914  fdm2d(       f, 0,   NC-1, 0 );
     915  fdm2d(     jac, 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 );
    916919}

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

    r1950 r1956  
    11SUBROUTINE calchim(ngrid,qy_c,declin,dtchim,            &
    2      ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
     2     ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc,zdyevapCH4)
    33
    44  !---------------------------------------------------------------------------------------------------------
    … …  
    5656  ! -----------------------------------------------------------------
    5757
     58  USE, INTRINSIC :: iso_c_binding
    5859  USE comchem_h
    5960  USE dimphy
    … …  
    8586
    8687  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT)  :: dqyc        ! Chemical species tendencies on GCM layers (mol/mol/s).
     88  REAL*8, DIMENSION(ngrid),           INTENT(OUT)  :: zdyevapCH4  ! Diagnostic surface methane pseudo-evaporation flux (mol/mol/s).
    8789
    8890  ! Local variables :
    … …  
    115117  REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
    116118  ! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
     119
     120  REAL(c_double) :: fluxCH4 ! Surface "evaporation" flux (mol/mol)
    117121
    118122  ! Saved variables initialized at firstcall
    … …  
    465469             nomqy_c,cqy,                               &
    466470             dtchim,latitude(ig)*180./pi,mass,md,       &
    467              kedd,botCH4,krate,reactif,                 &
     471             kedd,botCH4,fluxCH4,krate,reactif,         &
    468472             nom_prod,nom_perte,prod,perte,             &
    469473             aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,  &
    470474             htoh2,surfhaze)
     475
     476        zdyevapCH4(ig) = fluxCH4 / dtchim ! Diagnostic pseudo-evaporation ( due to readjustement to botCH4 value ) (mol/mol/s)
    471477
    472478        ! 5. Calculates tendencies on composition for advected tracers
    … …  
    488494
    489495
    490      ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
     496     ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
     497        zdyevapCH4(ig)  = zdyevapCH4(igm1)
    491498        dqyc(ig,:,:)    = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
    492499        ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations

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

    r1947 r1956  
    387387      real temp_eq(nlayer), press_eq(nlayer) ! Planetary averages for the init. of saturation profiles (K,mbar)
    388388
    389       ! Surface methane tank
    390       real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank (m)
    391 !$OMP THREADPRIVATE(tankCH4)
     389      ! Surface methane
     390      real, dimension(:), allocatable, save :: tankCH4    ! Depth of surface methane tank (m)
     391      real, dimension(:), allocatable, save :: zdyevapCH4 ! Surface pseudo-evaporation flux (chemistry keeping constant surface humidity) (mol/mol/s).
     392!$OMP THREADPRIVATE(tankCH4,zdyevapCH4)
    392393
    393394      ! -----******----- FOR MUPHYS OPTICS -----******-----
    … …  
    522523            allocate(dycchi(ngrid,nlayer,nkim)) ! only for chemical tracers
    523524            allocate(qysat(nlayer,nkim))
     525            allocate(zdyevapCH4(ngrid))
    524526           
    525527            ! Chemistry timestep
    … …  
    11991201                 ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
    12001202                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,  &
    1201                               zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
     1203                              zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi,zdyevapCH4)
    12021204               else ! 3D chemistry (or 1D run)
    12031205                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,  &
    1204                               pplay,pplev,zzlay,zzlev,dycchi)
     1206                              pplay,pplev,zzlay,zzlev,dycchi,zdyevapCH4)
    12051207               endif ! if moyzon
    12061208
    12071209            endif
     1210           
     1211            ! Add diagnostic-only surface pseudo-evapoaration in condensation tendency for bottom layer
     1212            zdqcond(:,1,chimi_indx(7)) = zdyevapCH4(:)*rat_mmol(chimi_indx(7))
    12081213           
    12091214            do iq=1,nkim
    … …  
    15591564             call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
    15601565           enddo
     1566           call writediagfi(ngrid,"fluxCH4","Surface CH4 pseudo-evaporation",'mol/mol/s',2,zdyevapCH4)
    15611567         endif
    15621568
Note: See TracChangeset for help on using the changeset viewer.