[2696] | 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 | |
---|
| 57 | return |
---|
| 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 |
---|