source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/exqssa.c @ 2697

Last change on this file since 2697 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 2.3 KB
Line 
1void INTEGRATE( double DT )
2{
3KPP_REAL P_VAR[NVAR], D_VAR[NVAR], V1[NVAR], V2[NVAR];
4int IsReject;
5KPP_REAL T, Tnext, STEP, STEPold, Told, SUP;
6KPP_REAL ERR, ERRold, ratio, factor, facmax, tmp;
7int i;
8
9  T = TIME;
10  Tnext = TIME + DT;
11  STEP = STEPMIN;
12  Told = T;
13  SUP  = 1e-14;
14  IsReject = 0;
15  ERR = .5;
16 
17/* -- BELOW THIS LIMIT USE TAYLOR INSTEAD OF EXP --- */
18
19  while ( T < Tnext ) {
20   
21    T = Told + STEP;
22    if ( T > Tnext ) {
23      STEP = Tnext - Told;
24      T = Tnext;
25    }
26
27    FSPLIT_VAR ( VAR,  P_VAR, D_VAR );
28
29    for( i = 0; i < NVAR; i++ ) {
30      if ( fabs(D_VAR[i]) > SUP ) {
31        ratio = P_VAR[i] / D_VAR[i];
32        tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) );
33        V1[i] = tmp * tmp * (VAR[i] - ratio) + ratio;
34        V2[i] = tmp * (VAR[i] - ratio) + ratio;
35      } else {
36        tmp = D_VAR[i] * STEP * 0.5;
37        V1[i] = VAR[i] + P_VAR[i] * STEP * ( 1 - tmp *
38               ( 1 - 2.0 / 3.0 * tmp ) );
39        V2[i] = VAR[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp *
40               ( 1 - 1.0 / 3.0 * tmp ) );
41      }
42    }
43
44    FSPLIT_VAR( V2,  P_VAR, D_VAR );
45
46    for( i = 0; i < NVAR; i++ ) {
47      if ( fabs(D_VAR[i]) > SUP ) {
48        ratio = P_VAR[i] / D_VAR[i];
49        tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) );
50        V2[i] = tmp * (V2[i] - ratio) + ratio;
51      } else {
52        tmp = D_VAR[i] * STEP * 0.5;
53        V2[i] = V2[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp *
54               ( 1 - 1.0 / 3.0 * tmp ) );
55      }
56    }
57/* ==== Extrapolation and error estimation ======== */
58
59    ERRold=ERR;
60    ERR=0.;
61    for( i = 0; i < NVAR; i++ ) {
62      V1[i] = 2.*V2[i] - V1[i];
63      tmp = (V2[i] - V1[i]) / (ATOL[i] + RTOL[i]*V2[i]); 
64      ERR = ERR + tmp*tmp;
65    }
66    ERR = sqrt(ERR/NVAR);
67    STEPold = STEP;
68
69/* ===== choosing the stepsize ==================== */
70
71    factor = 0.9 / pow(ERR,0.35) * pow(ERRold,0.2);
72    facmax = IsReject ? 1.0 : 5.0;
73
74    factor = max( 0.2, min(factor,facmax) );
75    STEP = min( STEPMAX, max(STEPMIN,factor*STEP) );
76
77/*================================================= */
78
79    if ( (ERR > 1) && (STEPold > STEPMIN) ) {
80      T = Told;
81      IsReject = 1;
82    } else {
83      IsReject = 0;
84      Told = T;
85      for( i = 0; i < NVAR; i++ ) 
86        VAR[i] = max( V1[i], 0.0 );
87      TIME = Tnext;
88    }
89  }
90}
Note: See TracBrowser for help on using the repository browser.