chem/KPP/kpp/kpp-2.1/drv/general.m

   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