detection/cougar_cycles.m

   1 function cycles(varargin)
   2 ! date
   3
   4 % base_datestr='2013-08-11 00:00:00';
   5 % for the cougar fire
   6 base_datestr='2015-08-11 00:00:00';
   7
   8 %multiplier to go to different time frame(s)
   9 mult = 0.5;
  10
  11 base=datenum(base_datestr);
  12 num_cycles=5;
  13 spinup_time=ones(1,num_cycles)*mult;
  14 cycle_length=ones(1,num_cycles)*mult;
  15 cycle_start =[0,2,4,6,8,10]*mult
  16 spinup_time =[2,2,2,2,2,2]*mult
  17 times_format='yyyy-mm-dd_HH:MM:SS';
  18 for i=1:num_cycles
  19     t(i).forecast_time=cycle_start(i+1)+1*mult;
  20     t(i).obs_start=cycle_start(i);
  21     t(i).obs_end=cycle_start(i+1)-1e-6;
  22     t(i).replay_start=cycle_start(i);
  23     t(i).replay_end=cycle_start(i+1);
  24     t(i).run_end=cycle_start(i+1)+3*mult;
  25     t(i).perimeter_time=t(i).replay_end*24*3600;
  26     forecast_times{i}=datestr(base+t(i).forecast_time,times_format);
  27     print_times(i)
  28 end
  29 print_times_table
  30
  31 if nargin,
  32     i=varargin{1};
  33     force=1;
  34 else
  35     i=input_num('cycle number',1);
  36     force=0;
  37 end
  38 %force = 1;
  39 link_namelist_command=sprintf('rm -f namelist.input; ln -s namelist.input_%i namelist.input',i);
  40 if i==0,
  41     fprintf('Initial simulation from %s to at least %s\n',base_datestr,forecast_times{1})
  42     fprintf('Execute %s now?\n',link_namelist_command);
  43     if input_num('0/1',1,force)
  44         if system(link_namelist_command),
  45              error('failed')
  46         end
  47     end
  48     disp('Run WRF-SFIRE and continue when done\n')
  49
  50 else
  51     print_times(i)
  52     system('ls -lh wrfout*')
  53     wrfout_time = base+t(i).forecast_time;
  54     wrfout{i}=['wrfout_d04_',datestr(wrfout_time,times_format)];
  55     if ~exist(wrfout{i},'file')
  56         fprintf('file %s does not exist\n',wrfout{i})
  57         wrfout_time = wrfout_time - 23.5/24;  % no wrfout produced on restart => written 30 min later
  58         wrfout{i}=['wrfout_d04_',datestr(wrfout_time,times_format)];
  59     end
  60     wrfrst{i}=['wrfrst_d04_',datestr(base+t(i).replay_start,times_format)];
  61     fprintf('%s %s %s %s\n','Reading fire arrival time at',forecast_times{i},' from ',wrfout{i})
  62     if t(i).replay_start==0;
  63        rewrite='wrfinput_d04';
  64        restart='.false.';
  65     else
  66        rewrite=wrfrst{i};
  67        restart='.true.';
  68     end
  69     fprintf('%s %s\n','Will write modified time into     ',rewrite)
  70     rewrite_bak=[rewrite,'.bak'];
  71     q=input_num(['1 to copy ',rewrite,' to ',rewrite_bak],1,force);
  72     if q,
  73        if system(['cp ',rewrite,' ',rewrite_bak]),
  74            warning('copy failed')
  75        end
  76     end
  77     w=read_wrfout_tign(wrfout{i},forecast_times{i});
  78     if isempty(w)
  79         display('Time step not found. Taking final timestep instead')
  80         w=read_wrfout_tign(wrfout{i}) % this will take final timestep
  81     end
  82     wrfout_bak=[wrfout{i},'.bak'];
  83     q=input_num(['1 to move ',wrfout{i},' to ',wrfout_bak],1,force);
  84     if q,
  85        movefile(wrfout{i},wrfout_bak);
  86     end
  87     time_bounds=[t(i).obs_start,t(i).obs_end,t(i).replay_start,t(i).replay_end]+base;
  88     savew=sprintf('w_%i',i);
  89     fprintf('saving to %s\n',savew)
  90     cycle=i;
  91     save(savew,'w','cycle','time_bounds','t')
  92     p=detect_fit_level2(cycle,time_bounds,[],w,force)
  93     fprintf('Saving p\n');
  94     pstr = sprintf('p_%i.mat',i);
  95     save(pstr,'p')
  96     print_times(i)
  97     fprintf('perimeter_time=%10.3f\nrestart=%s\n',t(i).perimeter_time,restart)
  98     q=sprintf('replace TIGN_G in %s and run\n %s\n [0/1]',rewrite,link_namelist_command);
  99     if input_num(q,1,force)
 100         ncreplace(rewrite,'TIGN_G',p.spinup)
 101         if system(link_namelist_command),
 102              error('link failed')
 103         end
 104     end
 105     disp('Run WRF-SFIRE and continue when done\n')
 106 end
 107
 108 function print_times(ii)
 109 ptime(ii,'Forecast used     ',t(ii).forecast_time)
 110 ptime(ii,'Observations start',t(ii).obs_start)
 111 ptime(ii,'Observations end  ',t(ii).obs_end)
 112 ptime(ii,'Replay start      ',t(ii).replay_start)
 113 ptime(ii,'Replay end        ',t(ii).replay_end)
 114 ptime(ii,'Run end           ',t(ii).run_end)
 115 fprintf('perimeter_time=%10.3f =  %d days \n',t(ii).perimeter_time,(t(ii).perimeter_time)/3600/24)
 116 end
 117
 118 function print_times_table
 119 fmt='  %g %g %g %g %g\n';
 120 fprintf(['Cycle             ',fmt],1:length(t))
 121 fprintf(['Forecast used     ',fmt],t(:).forecast_time)
 122 fprintf(['Observations start',fmt],t(:).obs_start)
 123 fprintf(['Observations end  ',fmt],t(:).obs_end)
 124 fprintf(['Replay start      ',fmt],t(:).replay_start)
 125 fprintf(['Replay end        ',fmt],t(:).replay_end)
 126 fprintf(['Run end           ',fmt],t(:).run_end)
 127 end
 128
 129 function ptime(ii,s,t)
 130         fprintf('Cycle %i %s%7.3f days %s\n',ii,s,t,datestr(t+base,'dd-mmm-yyyy HH:MM:SS'))
 131 end
 132
 133 end