source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general.m @ 4103

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

Merge of branch palm4u into trunk

File size: 1.9 KB
Line 
1
2  TSTART = 0;
3  TEND = TSTART + 600;
4  DT = 60.;
5  TEMP = 298;
6
7  RTOLS = 1.0e-6;
8  ATOLS = 1.0e-3;
9     
10  KPP_ROOT_Initialize;
11 
12  Options = odeset('AbsTol',ATOLS,'RelTol',RTOLS,'Jacobian',@KPP_ROOT_Jac_Chem);
13
14% ********** TIME LOOP *************************
15
16  C(1:KPP_NVAR) = VAR(1:KPP_NVAR);
17  C((KPP_NVAR+1):KPP_NSPEC) = FIX(1:KPP_NFIX);
18  DVAL = KPP_ROOT_GetMass( C );
19  if ( ~isempty(SMASS) )
20     fprintf('Initial Mass = %10.4e\n', DVAL(1:NMASS)/CFACTOR);
21  end   
22 
23%  KPP_ROOT_InitializeSaveData;
24
25%  disp(['Done[%] Time[h] ',SPC_NAMES(MONITOR(1:NMONITOR))])
26
27  TIME = TSTART;
28 
29  Tspan = linspace( TSTART, TEND, 100 );
30
31  [T, Y] = ode15s(@KPP_ROOT_Fun_Chem, Tspan, VAR, Options);
32
33  VAR(1:KPP_NVAR) = Y(length(T),1:KPP_NVAR)';
34  Y = [Y, ones(length(T),1)*FIX(:)'];
35 
36  C(1:KPP_NVAR) = VAR(1:KPP_NVAR);
37  C((NVAR+1):NSPEC) = FIX(1:NFIX);
38  DVAL = KPP_ROOT_GetMass( C );
39       
40  fprintf('done %6.1f,  %7.2h hours', (TIME-TSTART)/(TEND-TSTART), TIME/3600.);
41  disp( Y(:,MONITOR(1:NMONITOR))/CFACTOR );
42  if ( ~isempty(SMASS) )
43    fprintf('Final Mass = %10.4e\n', DVAL(1:NMASS)/CFACTOR);
44  end
45 
46  for i = 1:NMONITOR
47    figure; plot( (T)/3600, Y(:,MONITOR(i))/CFACTOR );
48    title( SPC_NAMES( MONITOR(i),:) ,'FontSize',12);
49    set(gca,'XLim',[TSTART,TEND]/3600,'FontSize',12);
50    xlabel('Time [ h ]','FontSize',12);
51    ylabel('Concentration','FontSize',12);
52  end 
53
54%   KPP_ROOT_FUNC_SaveData;
55%   KPP_ROOT_FUNC_CloseSaveData;
56
57return
58
59
60%  function P = KPP_ROOT_Fun_Chem(T, Y)
61%       
62%    global TIME FIX RCONST 
63%   
64%    Told = TIME;
65%    TIME = T;
66%    KPP_ROOT_Update_SUN;
67%    KPP_ROOT_Update_RCONST;
68%    P = KPP_ROOT_Fun( Y, FIX, RCONST );
69%    TIME = Told;
70%  return
71
72%   
73%  function J =  KPP_ROOT_Jac_Chem(T, Y)   
74%   
75%    global TIME FIX RCONST ;
76%   
77%    Told = TIME;
78%    TIME = T;
79%    KPP_ROOT_Update_SUN;
80%    KPP_ROOT_Update_RCONST;
81%    J = KPP_ROOT_Jac_SP( Y, FIX, RCONST );
82%    TIME = Told
83%  return             
Note: See TracBrowser for help on using the repository browser.