standalone/run_moisture_test.m

   1 function run_moisture_test(varargin)
   2 if nargin>=1,
   3     ncfile=varargin{1};
   4 else
   5     ncfile=[];
   6 end
   7 if isempty(ncfile),
   8     ncfile='moisture100.nc';
   9 end
  10 if nargin>=2,
  11     hr=varargin{2};
  12 else
  13     hr=[0,realmax];
  14 end
  15 ncload(ncfile)
  16 Times=ncread(ncfile,'Times');
  17 [d,hours]=make_moisture_input(Times,t2,psfc,q2,rain);
  18 setenv('DYLD_LIBRARY_PATH','/opt/local/lib:/opt/local/lib/gcc46')
  19 system('./moisture_test.exe');
  20 load moisture_output.txt
  21 out.t2   =moisture_output(:,3);
  22 out.t2_mid=midpoints(out.t2);
  23 out.psfc =moisture_output(:,4);
  24 out.q2   =moisture_output(:,5);
  25 % rh_fire is computed from midpoint values in the fortran code
  26 out.rain=moisture_output(:,6);
  27 out.rh_fire=moisture_output(:,7);
  28 out.fmc_equi =moisture_output(:,8);
  29 out.tlag     =moisture_output(:,9);
  30 out.fmc_gc   =moisture_output(:,10);
  31 out.id='from moisture_test.exe';
  32 out.id='';
  33 out.hours=hours;
  34 out.hours_mid=midpoints(hours);
  35 out.r=find(hours>=hr(1) & hours<=hr(2)); % range to display
  36 orig.t2=t2(~d);
  37 orig.t2_mid=midpoints(orig.t2);
  38 orig.psfc=psfc(~d);
  39 orig.q2=q2(~d);
  40 orig.rain=rain(~d);
  41 orig.rh_fire=rh_fire(~d);
  42 orig.fmc_equi=fmc_equi(1,~d)';
  43 orig.fmc_gc=fmc_gc(1,~d)';
  44 orig.id=['from ',ncfile];
  45 orig.hours=hours;
  46 orig.hours_mid=midpoints(hours);
  47 orig.r=out.r;
  48 figure(1)
  49 plot(hours,orig.fmc_gc,'r--',hours,out.fmc_gc,'k-')
  50 xlabel hours
  51 ylabel kg/kg
  52 title('Fuel moisture')
  53 h=legend(orig.id,out.id);
  54 set(h,'Interpreter','none')
  55 figure(2)
  56 plot(hours,orig.fmc_equi,'r--',hours,out.fmc_equi,'k-')
  57 xlabel hours
  58 set(ylabel('kg/kg'),'Interpreter','none')
  59 title('Equilibrium moisture')
  60 h=legend(orig.id,out.id);
  61 set(h,'Interpreter','none')
  62 plot_moisture(3,orig)
  63 plot_moisture(4,out)
  64 plot_all(5,orig)
  65 plot_all(6,out)
  66 end
  67
  68 function plot_moisture(f,s)
  69 figure(f);
  70 [d,w]=equilibrium_moisture(s.rh_fire,s.t2_mid);
  71 plot(s.hours_mid(s.r),d(s.r),'b--',...
  72      s.hours_mid(s.r),w(s.r),'b-.',...
  73      s.hours_mid(s.r),s.fmc_equi(s.r),'r--',...
  74      s.hours(s.r),s.fmc_gc(s.r),'k-')
  75 xlabel hours
  76 set(ylabel('Fuel moisture contents (kg / kg)'),'Interpreter','none')
  77 h=title(['Fuel moisture ',s.id]);
  78 set(h,'Interpreter','none')
  79 h=legend('Drying','Wetting','Equilibrium','Actual');
  80 set(h,'Interpreter','none')
  81 setvmax(0.3)
  82 end
  83
  84 function plot_all(f,s)
  85 figure(f)
  86 [d,w]=equilibrium_moisture(s.rh_fire,s.t2_mid);
  87 plot(s.hours(s.r),s.t2(s.r)/max(s.t2),...
  88     s.hours(s.r),s.psfc(s.r)/max(s.psfc),...
  89     s.hours(s.r),s.q2(s.r)/max(s.q2),...
  90     s.hours(s.r),s.rh_fire(s.r),...
  91     s.hours_mid(s.r),s.fmc_equi(s.r),...
  92     s.hours_mid(s.r),d(s.r),'g-.',...
  93     s.hours_mid(s.r),w(s.r),'b-.',...
  94     s.hours(s.r),s.fmc_gc(s.r))
  95 h=title(['All variables ',s.id]);
  96 set(h,'Interpreter','none');
  97 xlabel hours
  98 h=legend('Scaled T2',...
  99     'Scaled PSFC',...
 100     'Scaled Q2',...
 101     'RH',...
 102     'FMC_EQUI',...
 103     'Drying equilibrium',...
 104     'Wetting equilibrium',...
 105     'FMC_GC')
 106 set(h,'Interpreter','none');
 107 setvmax(1)
 108 end
 109
 110 function setvmax(v)
 111 a=axis;
 112 a(3:4)=[0,v];
 113 axis(a)
 114 end
 115
 116 function T=midpoints(T2)
 117 T=[T2(1);0.5*(T2(1:end-1)+T2(2:end))];
 118 end