- 积分
- 38
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-1-2
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2015-1-14 17:10:54
|
显示全部楼层
% clear everything
clear all;
clc;
path(path,'D:\DataProgram\MatlabLib');
path(path,'D:\DataProgram\MatlabLib\savefigure');
% global 1 by 1 degree data
lon = -179.5:1.0:179.5;
lat = -89.5:1.0:89.5;
n_lon = 360;
n_lat = 180;
% find lon and lat index for amazon, africa, southeast asia, southeast
% america, europe
% (65.5W, 7.5S) Amazon
%( 23.5E,2.5S) Africa,
% (113.5E,1.5N) Southeast Asia
% (83.5W, 34.5N) Southeast America
% (4.5E, 45.5N) in Europe
lon_amazon = -65.5; lat_amazon = -7.5;
lon_amazon_index = lon == lon_amazon;
lat_amazon_index = lat == lat_amazon;
lon_africa = 23.5; lat_africa = -2.5;
lon_africa_index = lon == lon_africa;
lat_africa_index = lat == lat_africa;
lon_asia = 113.5; lat_asia = 1.5;
lon_asia_index = lon == lon_asia;
lat_asia_index = lat == lat_asia;
lon_america = -83.5; lat_america = 34.5;
lon_america_index = lon == lon_america;
lat_america_index = lat == lat_america;
lon_europe = 4.5; lat_europe = 45.5;
lon_europe_index = lon == lon_europe;
lat_europe_index = lat == lat_europe;
% define input and output directories
inputdir = 'D:\Data\CLM_output\Monthly';
outputdir = 'D:\DataProgram\CLM_output';
% check if input directories exists
if (exist(inputdir,'dir') == 0)
error('ErrorMsg:Inputdir', ...
'Input directory %s does not exist.', inputdir);
end;
% check if input directories exists. if not, make one
if (exist(outputdir,'dir') == 0)
str = sprintf('Output directory %s does not exist.', ...
outputdir);
disp(str);
disp('The program is making for you.');
[istatus,mess,messid] = mkdir(outputdir);
if (istatus == 1)
disp('The program has made the directory sucessfully');
else
error('ErrorMsg:MakeOutputdir', ...
'The program can not make the directory %s.', outputdir);
end
end
% ======================================================
% ctl: the CLM3 default scheme
% ======================================================
inputfiles = dir([inputdir, '\', 'CLM3_CTL_60min.clm2.h0.2002*.nc']);
if (length(inputfiles) < 1)
error('ErrorMsg:InputdirFile', ...
'There is no required file in the input directory %s.', inputdir);
end
depth = [18 28 45 75 124 204 336 554 913 1137];
% seasonal
for i = 1:length(inputfiles)
indirfile = [inputdir, '\', inputfiles(i).name];
if (exist(indirfile,'file') == 2)
ncid = netcdf.open(indirfile,'NC_NOWRITE');
if ncid > 0
% rainfall (mm/s) --> (mm/day)
rainid = netcdf.inqVarID(ncid,'RAIN');
rain = netcdf.getVar(ncid,rainid);
rain = double(24.0*3600.0*rain);
rain(rain > 2.0e+30) = NaN;
rain_amazon(i) = rain(lon_amazon_index,lat_amazon_index);
rain_africa(i) = rain(lon_africa_index,lat_africa_index);
rain_asia(i) = rain(lon_asia_index,lat_asia_index);
rain_america(i) = rain(lon_america_index,lat_america_index);
rain_europe(i) = rain(lon_europe_index,lat_europe_index);
% incident solar radiation W/m2
fsdsid = netcdf.inqVarID(ncid,'FSDS');
fsds = netcdf.getVar(ncid,fsdsid);
fsds = double(24.0*3600.0*fsds);
fsds(fsds > 2.0e+30) = NaN;
fsds_amazon(i) = fsds(lon_amazon_index,lat_amazon_index);
fsds_africa(i) = fsds(lon_africa_index,lat_africa_index);
fsds_asia(i) = fsds(lon_asia_index,lat_asia_index);
fsds_america(i) = fsds(lon_america_index,lat_america_index);
fsds_europe(i) = fsds(lon_europe_index,lat_europe_index);
% LAI
elaiid = netcdf.inqVarID(ncid,'ELAI');
elai = netcdf.getVar(ncid,elaiid);
elai = double(elai);
elai(elai > 2.0e+30) = NaN;
elai_amazon(i) = elai(lon_amazon_index,lat_amazon_index);
elai_africa(i) = elai(lon_africa_index,lat_africa_index);
elai_asia(i) = elai(lon_asia_index,lat_asia_index);
elai_america(i) = elai(lon_america_index,lat_america_index);
elai_europe(i) = elai(lon_europe_index,lat_europe_index);
% Sensible heat flux (W/m2)
fshid = netcdf.inqVarID(ncid,'FSH');
fsh = netcdf.getVar(ncid,fshid);
fsh = double(fsh);
fsh(fsh > 2.0e+30) = NaN;
fsh_amazon_ctl(i) = fsh(lon_amazon_index,lat_amazon_index);
fsh_africa_ctl(i) = fsh(lon_africa_index,lat_africa_index);
fsh_asia_ctl(i) = fsh(lon_asia_index,lat_asia_index);
fsh_america_ctl(i) = fsh(lon_america_index,lat_america_index);
fsh_europe_ctl(i) = fsh(lon_europe_index,lat_europe_index);
% heat flux for plant transpiration (W/m2)
fctrid = netcdf.inqVarID(ncid,'FCTR');
fctr = netcdf.getVar(ncid,fctrid);
fctr = double(fctr);
fctr(fctr > 2.0e+30) = NaN;
% heat flux for interception loss (W/m2)
fcevid = netcdf.inqVarID(ncid,'FCEV');
fcev = netcdf.getVar(ncid,fcevid);
fcev = double(fcev);
fcev(fcev > 2.0e+30) = NaN;
% heat flux for soil evaporation (W/m2)
fgevid = netcdf.inqVarID(ncid,'FGEV');
fgev = netcdf.getVar(ncid,fgevid);
fgev = double(fgev);
fgev(fgev > 2.0e+30) = NaN;
% latent heat flux (W/m2)
fle = fctr + fcev + fgev;
fle_amazon_ctl(i) = fle(lon_amazon_index,lat_amazon_index);
fle_africa_ctl(i) = fle(lon_africa_index,lat_africa_index);
fle_asia_ctl(i) = fle(lon_asia_index,lat_asia_index);
fle_america_ctl(i) = fle(lon_america_index,lat_america_index);
fle_europe_ctl(i) = fle(lon_europe_index,lat_europe_index);
% canopy interception (mm/s) --> (mm/day)
qvegeid = netcdf.inqVarID(ncid,'QVEGE');
qvege = netcdf.getVar(ncid,qvegeid);
qvege = double(24.0*3600.0*qvege);
qvege(qvege > 2.0e+30) = NaN;
qvege_amazon_ctl(i) = qvege(lon_amazon_index,lat_amazon_index);
qvege_africa_ctl(i) = qvege(lon_africa_index,lat_africa_index);
qvege_asia_ctl(i) = qvege(lon_asia_index,lat_asia_index);
qvege_america_ctl(i) = qvege(lon_america_index,lat_america_index);
qvege_europe_ctl(i) = qvege(lon_europe_index,lat_europe_index);
% plant transpiration (mm/s) --> (mm/day)
qvegtid = netcdf.inqVarID(ncid,'QVEGT');
qvegt = netcdf.getVar(ncid,qvegtid);
qvegt = double(24.0*3600.0*qvegt);
qvegt(qvegt > 2.0e+30) = NaN;
qvegt_amazon_ctl(i) = qvegt(lon_amazon_index,lat_amazon_index);
qvegt_africa_ctl(i) = qvegt(lon_africa_index,lat_africa_index);
qvegt_asia_ctl(i) = qvegt(lon_asia_index,lat_asia_index);
qvegt_america_ctl(i) = qvegt(lon_america_index,lat_america_index);
qvegt_europe_ctl(i) = qvegt(lon_europe_index,lat_europe_index);
% soil evaporation (mm/s) --> (mm/day)
qsoilid = netcdf.inqVarID(ncid,'QSOIL');
qsoil = netcdf.getVar(ncid,qsoilid);
qsoil = double(24.0*3600.0*qsoil);
qsoil(qsoil > 2.0e+30) = NaN;
qsoil_amazon_ctl(i) = qsoil(lon_amazon_index,lat_amazon_index);
qsoil_africa_ctl(i) = qsoil(lon_africa_index,lat_africa_index);
qsoil_asia_ctl(i) = qsoil(lon_asia_index,lat_asia_index);
qsoil_america_ctl(i) = qsoil(lon_america_index,lat_america_index);
qsoil_europe_ctl(i) = qsoil(lon_europe_index,lat_europe_index);
% total evapotranspiration (mm/s) --> (mm/day)
et = qvege + qvegt + qsoil;
et_amazon_ctl(i) = et(lon_amazon_index,lat_amazon_index);
et_africa_ctl(i) = et(lon_africa_index,lat_africa_index);
et_asia_ctl(i) = et(lon_asia_index,lat_asia_index);
et_america_ctl(i) = et(lon_america_index,lat_america_index);
et_europe_ctl(i) = et(lon_europe_index,lat_europe_index);
% surface runoff (mm/s) --> (mm/day)
qoverid = netcdf.inqVarID(ncid,'QOVER');
qover = netcdf.getVar(ncid,qoverid);
qover = double(24.0*3600.0*qover);
qover(qover > 2.0e+30) = NaN;
qover_amazon_ctl(i) = qover(lon_amazon_index,lat_amazon_index);
qover_africa_ctl(i) = qover(lon_africa_index,lat_africa_index);
qover_asia_ctl(i) = qover(lon_asia_index,lat_asia_index);
qover_america_ctl(i) = qover(lon_america_index,lat_america_index);
qover_europe_ctl(i) = qover(lon_europe_index,lat_europe_index);
% sub-surface runoff (mm/s) --> (mm/day)
qdraiid = netcdf.inqVarID(ncid,'QDRAI');
qdrai = netcdf.getVar(ncid,qdraiid);
qdrai = double(24.0*3600.0*qdrai);
qdrai(qdrai > 2.0e+30) = NaN;
qdrai_amazon_ctl(i) = qdrai(lon_amazon_index,lat_amazon_index);
qdrai_africa_ctl(i) = qdrai(lon_africa_index,lat_africa_index);
qdrai_asia_ctl(i) = qdrai(lon_asia_index,lat_asia_index);
qdrai_america_ctl(i) = qdrai(lon_america_index,lat_america_index);
qdrai_europe_ctl(i) = qdrai(lon_europe_index,lat_europe_index);
% runoff (mm/s) --> (mm/day)
runoff = qover + qdrai;
runoff_amazon_ctl(i) = runoff(lon_amazon_index,lat_amazon_index);
runoff_africa_ctl(i) = runoff(lon_africa_index,lat_africa_index);
runoff_asia_ctl(i) = runoff(lon_asia_index,lat_asia_index);
runoff_america_ctl(i) = runoff(lon_america_index,lat_america_index);
runoff_europe_ctl(i) = runoff(lon_europe_index,lat_europe_index);
% ground temperature (K)
tgid = netcdf.inqVarID(ncid,'TG');
tg = netcdf.getVar(ncid,tgid);
tg = double(tg);
tg(tg > 2.0e+30) = NaN;
tg_amazon_ctl(i) = tg(lon_amazon_index,lat_amazon_index);
tg_africa_ctl(i) = tg(lon_africa_index,lat_africa_index);
tg_asia_ctl(i) = tg(lon_asia_index,lat_asia_index);
tg_america_ctl(i) = tg(lon_america_index,lat_america_index);
tg_europe_ctl(i) = tg(lon_europe_index,lat_europe_index);
% vegetation temperature (K)
tvid = netcdf.inqVarID(ncid,'TV');
tv = netcdf.getVar(ncid,tvid);
tv = double(tv);
tv(tv > 2.0e+30) = NaN;
tv_amazon_ctl(i) = tv(lon_amazon_index,lat_amazon_index);
tv_africa_ctl(i) = tv(lon_africa_index,lat_africa_index);
tv_asia_ctl(i) = tv(lon_asia_index,lat_asia_index);
tv_america_ctl(i) = tv(lon_america_index,lat_america_index);
tv_europe_ctl(i) = tv(lon_europe_index,lat_europe_index);
% soil water (kg/m2) == mm
soilliqid = netcdf.inqVarID(ncid,'SOILLIQ');
soilliq = netcdf.getVar(ncid,soilliqid);
soilliq = double(soilliq);
soilliq(soilliq > 2.0e+30) = NaN;
soiliceid = netcdf.inqVarID(ncid,'SOILICE');
soilice = netcdf.getVar(ncid,soiliceid);
soilice = double(soilice);
soilice(soilice > 2.0e+30) = NaN;
soil = soilliq + soilice;
soil(soil > 2.0e+30) = NaN;
% total soil moisture for 10 layers (3.4m)
soil_total = sum(soil,3,'native');
soil_amazon_ctl(i) = soil_total(lon_amazon_index,lat_amazon_index);
soil_africa_ctl(i) = soil_total(lon_africa_index,lat_africa_index);
soil_asia_ctl(i) = soil_total(lon_asia_index,lat_asia_index);
soil_america_ctl(i) = soil_total(lon_america_index,lat_america_index);
soil_europe_ctl(i) = soil_total(lon_europe_index,lat_europe_index);
soil_1m_total = sum(soil(:,:,1:7),3,'native') + soil(:,:,8)*0.1711;
soil_1m_amazon_ctl(i) = soil_1m_total(lon_amazon_index,lat_amazon_index);
soil_1m_africa_ctl(i) = soil_1m_total(lon_africa_index,lat_africa_index);
soil_1m_asia_ctl(i) = soil_1m_total(lon_asia_index,lat_asia_index);
soil_1m_america_ctl(i) = soil_1m_total(lon_america_index,lat_america_index);
soil_1m_europe_ctl(i) = soil_1m_total(lon_europe_index,lat_europe_index);
% soil temperature (K)
tsoiid = netcdf.inqVarID(ncid,'TSOI');
tsoi = netcdf.getVar(ncid,tsoiid);
tsoi = double(tsoi);
tsoi(tsoi > 2.0e+30) = NaN;
tsoi_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
+ tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
+ tsoi(:,:,7)*336 + tsoi(:,:,8)*554 + tsoi(:,:,9)*913 + tsoi(:,:,9)*1137)/3400;
tsoi_amazon_ctl(i) = tsoi_ave(lon_amazon_index,lat_amazon_index);
tsoi_africa_ctl(i) = tsoi_ave(lon_africa_index,lat_africa_index);
tsoi_asia_ctl(i) = tsoi_ave(lon_asia_index,lat_asia_index);
tsoi_america_ctl(i) = tsoi_ave(lon_america_index,lat_america_index);
tsoi_europe_ctl(i) = tsoi_ave(lon_europe_index,lat_europe_index);
tsoi_1m_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
+ tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
+ tsoi(:,:,7)*336 + tsoi(:,:,8)*171)/1000;
tsoi_1m_amazon_ctl(i) = tsoi_1m_ave(lon_amazon_index,lat_amazon_index);
tsoi_1m_africa_ctl(i) = tsoi_1m_ave(lon_africa_index,lat_africa_index);
tsoi_1m_asia_ctl(i) = tsoi_1m_ave(lon_asia_index,lat_asia_index);
tsoi_1m_america_ctl(i) = tsoi_1m_ave(lon_america_index,lat_america_index);
tsoi_1m_europe_ctl(i) = tsoi_1m_ave(lon_europe_index,lat_europe_index);
netcdf.close(ncid);
else
error('ErrorMsg:indirfileOpen', ...
'File %s can not be opened.', indirfile);
end
else
error('ErrorMsg:indirfileExistOrNot', ...
'File %s does not exist.', indirfile);
end
end
% ======================================================
% v2: the improved scheme
% ======================================================
inputfiles = dir([inputdir, '\', 'CLM3_EXP3_v2_60min.clm2.h0.2002*.nc']);
if (length(inputfiles) < 1)
error('ErrorMsg:InputdirFile', ...
'There is no required file in the input directory %s.', inputdir);
end
% seasonal
for i = 1:length(inputfiles)
indirfile = [inputdir, '\', inputfiles(i).name];
if (exist(indirfile,'file') == 2)
ncid = netcdf.open(indirfile,'NC_NOWRITE');
if ncid > 0
% rainfall (mm/s) --> (mm/day)
% rainid = netcdf.inqVarID(ncid,'RAIN');
% rain = netcdf.getVar(ncid,rainid);
% rain = double(24.0*3600.0*rain);
% rain(rain > 2.0e+30) = NaN;
% rain_amazon(i) = rain(lon_amazon_index,lat_amazon_index);
% rain_africa(i) = rain(lon_africa_index,lat_africa_index);
% rain_asia(i) = rain(lon_asia_index,lat_asia_index);
% rain_america(i) = rain(lon_america_index,lat_america_index);
% rain_europe(i) = rain(lon_europe_index,lat_europe_index);
%
% % LAI
% elaiid = netcdf.inqVarID(ncid,'ELAI');
% elai = netcdf.getVar(ncid,elaiid);
% elai = double(elai);
% elai(elai > 2.0e+30) = NaN;
% elai_amazon(i) = elai(lon_amazon_index,lat_amazon_index);
% elai_africa(i) = elai(lon_africa_index,lat_africa_index);
% elai_asia(i) = elai(lon_asia_index,lat_asia_index);
% elai_america(i) = elai(lon_america_index,lat_america_index);
% elai_europe(i) = elai(lon_europe_index,lat_europe_index);
% Sensible heat flux (W/m2)
fshid = netcdf.inqVarID(ncid,'FSH');
fsh = netcdf.getVar(ncid,fshid);
fsh = double(fsh);
fsh(fsh > 2.0e+30) = NaN;
fsh_amazon_exp(i) = fsh(lon_amazon_index,lat_amazon_index);
fsh_africa_exp(i) = fsh(lon_africa_index,lat_africa_index);
fsh_asia_exp(i) = fsh(lon_asia_index,lat_asia_index);
fsh_america_exp(i) = fsh(lon_america_index,lat_america_index);
fsh_europe_exp(i) = fsh(lon_europe_index,lat_europe_index);
% heat flux for plant transpiration (W/m2)
fctrid = netcdf.inqVarID(ncid,'FCTR');
fctr = netcdf.getVar(ncid,fctrid);
fctr = double(fctr);
fctr(fctr > 2.0e+30) = NaN;
% heat flux for interception loss (W/m2)
fcevid = netcdf.inqVarID(ncid,'FCEV');
fcev = netcdf.getVar(ncid,fcevid);
fcev = double(fcev);
fcev(fcev > 2.0e+30) = NaN;
% heat flux for soil evaporation (W/m2)
fgevid = netcdf.inqVarID(ncid,'FGEV');
fgev = netcdf.getVar(ncid,fgevid);
fgev = double(fgev);
fgev(fgev > 2.0e+30) = NaN;
% latent heat flux (W/m2)
fle = fctr + fcev + fgev;
fle_amazon_exp(i) = fle(lon_amazon_index,lat_amazon_index);
fle_africa_exp(i) = fle(lon_africa_index,lat_africa_index);
fle_asia_exp(i) = fle(lon_asia_index,lat_asia_index);
fle_america_exp(i) = fle(lon_america_index,lat_america_index);
fle_europe_exp(i) = fle(lon_europe_index,lat_europe_index);
% canopy interception (mm/s) --> (mm/day)
qvegeid = netcdf.inqVarID(ncid,'QVEGE');
qvege = netcdf.getVar(ncid,qvegeid);
qvege = double(24.0*3600.0*qvege);
qvege(qvege > 2.0e+30) = NaN;
qvege_amazon_exp(i) = qvege(lon_amazon_index,lat_amazon_index);
qvege_africa_exp(i) = qvege(lon_africa_index,lat_africa_index);
qvege_asia_exp(i) = qvege(lon_asia_index,lat_asia_index);
qvege_america_exp(i) = qvege(lon_america_index,lat_america_index);
qvege_europe_exp(i) = qvege(lon_europe_index,lat_europe_index);
% plant transpiration (mm/s) --> (mm/day)
qvegtid = netcdf.inqVarID(ncid,'QVEGT');
qvegt = netcdf.getVar(ncid,qvegtid);
qvegt = double(24.0*3600.0*qvegt);
qvegt(qvegt > 2.0e+30) = NaN;
qvegt_amazon_exp(i) = qvegt(lon_amazon_index,lat_amazon_index);
qvegt_africa_exp(i) = qvegt(lon_africa_index,lat_africa_index);
qvegt_asia_exp(i) = qvegt(lon_asia_index,lat_asia_index);
qvegt_america_exp(i) = qvegt(lon_america_index,lat_america_index);
qvegt_europe_exp(i) = qvegt(lon_europe_index,lat_europe_index);
% soil evaporation (mm/s) --> (mm/day)
qsoilid = netcdf.inqVarID(ncid,'QSOIL');
qsoil = netcdf.getVar(ncid,qsoilid);
qsoil = double(24.0*3600.0*qsoil);
qsoil(qsoil > 2.0e+30) = NaN;
qsoil_amazon_exp(i) = qsoil(lon_amazon_index,lat_amazon_index);
qsoil_africa_exp(i) = qsoil(lon_africa_index,lat_africa_index);
qsoil_asia_exp(i) = qsoil(lon_asia_index,lat_asia_index);
qsoil_america_exp(i) = qsoil(lon_america_index,lat_america_index);
qsoil_europe_exp(i) = qsoil(lon_europe_index,lat_europe_index);
% total evapotranspiration (mm/s) --> (mm/day)
et = qvege + qvegt + qsoil;
et_amazon_exp(i) = et(lon_amazon_index,lat_amazon_index);
et_africa_exp(i) = et(lon_africa_index,lat_africa_index);
et_asia_exp(i) = et(lon_asia_index,lat_asia_index);
et_america_exp(i) = et(lon_america_index,lat_america_index);
et_europe_exp(i) = et(lon_europe_index,lat_europe_index);
% surface runoff (mm/s) --> (mm/day)
qoverid = netcdf.inqVarID(ncid,'QOVER');
qover = netcdf.getVar(ncid,qoverid);
qover = double(24.0*3600.0*qover);
qover(qover > 2.0e+30) = NaN;
qover_amazon_exp(i) = qover(lon_amazon_index,lat_amazon_index);
qover_africa_exp(i) = qover(lon_africa_index,lat_africa_index);
qover_asia_exp(i) = qover(lon_asia_index,lat_asia_index);
qover_america_exp(i) = qover(lon_america_index,lat_america_index);
qover_europe_exp(i) = qover(lon_europe_index,lat_europe_index);
% sub-surface runoff (mm/s) --> (mm/day)
qdraiid = netcdf.inqVarID(ncid,'QDRAI');
qdrai = netcdf.getVar(ncid,qdraiid);
qdrai = double(24.0*3600.0*qdrai);
qdrai(qdrai > 2.0e+30) = NaN;
qdrai_amazon_exp(i) = qdrai(lon_amazon_index,lat_amazon_index);
qdrai_africa_exp(i) = qdrai(lon_africa_index,lat_africa_index);
qdrai_asia_exp(i) = qdrai(lon_asia_index,lat_asia_index);
qdrai_america_exp(i) = qdrai(lon_america_index,lat_america_index);
qdrai_europe_exp(i) = qdrai(lon_europe_index,lat_europe_index);
% runoff (mm/s) --> (mm/day)
runoff = qover + qdrai;
runoff_amazon_exp(i) = runoff(lon_amazon_index,lat_amazon_index);
runoff_africa_exp(i) = runoff(lon_africa_index,lat_africa_index);
runoff_asia_exp(i) = runoff(lon_asia_index,lat_asia_index);
runoff_america_exp(i) = runoff(lon_america_index,lat_america_index);
runoff_europe_exp(i) = runoff(lon_europe_index,lat_europe_index);
% ground temperature (K)
tgid = netcdf.inqVarID(ncid,'TG');
tg = netcdf.getVar(ncid,tgid);
tg = double(tg);
tg(tg > 2.0e+30) = NaN;
tg_amazon_exp(i) = tg(lon_amazon_index,lat_amazon_index);
tg_africa_exp(i) = tg(lon_africa_index,lat_africa_index);
tg_asia_exp(i) = tg(lon_asia_index,lat_asia_index);
tg_america_exp(i) = tg(lon_america_index,lat_america_index);
tg_europe_exp(i) = tg(lon_europe_index,lat_europe_index);
% vegetation temperature (K)
tvid = netcdf.inqVarID(ncid,'TV');
tv = netcdf.getVar(ncid,tvid);
tv = double(tv);
tv(tv > 2.0e+30) = NaN;
tv_amazon_exp(i) = tv(lon_amazon_index,lat_amazon_index);
tv_africa_exp(i) = tv(lon_africa_index,lat_africa_index);
tv_asia_exp(i) = tv(lon_asia_index,lat_asia_index);
tv_america_exp(i) = tv(lon_america_index,lat_america_index);
tv_europe_exp(i) = tv(lon_europe_index,lat_europe_index);
% soil water (kg/m2) == mm
soilliqid = netcdf.inqVarID(ncid,'SOILLIQ');
soilliq = netcdf.getVar(ncid,soilliqid);
soilliq = double(soilliq);
soilliq(soilliq > 2.0e+30) = NaN;
soiliceid = netcdf.inqVarID(ncid,'SOILICE');
soilice = netcdf.getVar(ncid,soiliceid);
soilice = double(soilice);
soilice(soilice > 2.0e+30) = NaN;
soil = soilliq + soilice;
soil(soil > 2.0e+30) = NaN;
% total soil moisture for 10 layers (3.4m)
soil_total = sum(soil,3,'native');
soil_amazon_exp(i) = soil_total(lon_amazon_index,lat_amazon_index);
soil_africa_exp(i) = soil_total(lon_africa_index,lat_africa_index);
soil_asia_exp(i) = soil_total(lon_asia_index,lat_asia_index);
soil_america_exp(i) = soil_total(lon_america_index,lat_america_index);
soil_europe_exp(i) = soil_total(lon_europe_index,lat_europe_index);
soil_1m_total = sum(soil(:,:,1:7),3,'native') + soil(:,:,8)*0.1711;
soil_1m_amazon_exp(i) = soil_1m_total(lon_amazon_index,lat_amazon_index);
soil_1m_africa_exp(i) = soil_1m_total(lon_africa_index,lat_africa_index);
soil_1m_asia_exp(i) = soil_1m_total(lon_asia_index,lat_asia_index);
soil_1m_america_exp(i) = soil_1m_total(lon_america_index,lat_america_index);
soil_1m_europe_exp(i) = soil_1m_total(lon_europe_index,lat_europe_index);
% soil temperature (K)
tsoiid = netcdf.inqVarID(ncid,'TSOI');
tsoi = netcdf.getVar(ncid,tsoiid);
tsoi = double(tsoi);
tsoi(tsoi > 2.0e+30) = NaN;
tsoi_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
+ tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
+ tsoi(:,:,7)*336 + tsoi(:,:,8)*554 + tsoi(:,:,9)*913 + tsoi(:,:,9)*1137)/3400;
tsoi_amazon_exp(i) = tsoi_ave(lon_amazon_index,lat_amazon_index);
tsoi_africa_exp(i) = tsoi_ave(lon_africa_index,lat_africa_index);
tsoi_asia_exp(i) = tsoi_ave(lon_asia_index,lat_asia_index);
tsoi_america_exp(i) = tsoi_ave(lon_america_index,lat_america_index);
tsoi_europe_exp(i) = tsoi_ave(lon_europe_index,lat_europe_index);
tsoi_1m_ave = (tsoi(:,:,1)*18 + tsoi(:,:,2)*28 + tsoi(:,:,3)*45 ...
+ tsoi(:,:,4)*75 + tsoi(:,:,5)*124 + tsoi(:,:,6)*204 ...
+ tsoi(:,:,7)*336 + tsoi(:,:,8)*171)/1000;
tsoi_1m_amazon_exp(i) = tsoi_1m_ave(lon_amazon_index,lat_amazon_index);
tsoi_1m_africa_exp(i) = tsoi_1m_ave(lon_africa_index,lat_africa_index);
tsoi_1m_asia_exp(i) = tsoi_1m_ave(lon_asia_index,lat_asia_index);
tsoi_1m_america_exp(i) = tsoi_1m_ave(lon_america_index,lat_america_index);
tsoi_1m_europe_exp(i) = tsoi_1m_ave(lon_europe_index,lat_europe_index);
netcdf.close(ncid);
else
error('ErrorMsg:indirfileOpen', ...
'File %s can not be opened.', indirfile);
end
else
error('ErrorMsg:indirfileExistOrNot', ...
'File %s does not exist.', indirfile);
end
end
% rain coverage
inputdir = 'D:\Data\CLM_output\ForcingAndObservedRain';
inputfiles = dir([inputdir, '\', 'raincoverage-*.txt']);
if (length(inputfiles) < 1)
error('ErrorMsg:InputdirFile', ...
'There is no required file in the input directory %s.', inputdir);
end
for i = 1:length(inputfiles)
indirfile = [inputdir, '\', inputfiles(i).name];
fid = fopen(indirfile,'r');
raincoverage = fscanf(fid, '%f', [360, 180]);
raincoverage_amazon(i) = raincoverage(lon_amazon_index,lat_amazon_index);
raincoverage_africa(i) = raincoverage(lon_africa_index,lat_africa_index);
raincoverage_asia(i) = raincoverage(lon_asia_index,lat_asia_index);
raincoverage_america(i) = raincoverage(lon_america_index,lat_america_index);
raincoverage_europe(i) = raincoverage(lon_europe_index,lat_europe_index);
fclose(fid);
end
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% &&&&&&&&&&&&&&&& plot &&&&&&&&&&&&&&&&&&&&&&&&&&&&
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
xx = 0:1:11;
scrsz = get(0,'ScreenSize');
figure('Position',[5 55 scrsz(3)-10 scrsz(4)-115]);
subplot1(2,5,'Gap',[0.04 0.07],'XTickL','Margin','YTickL','All','FontS',8,'Min',[0.03 0.1],'Max',[0.99 0.99]);
subplot1(1);
plot(xx,rain_amazon,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
ylabel('Rainfall (mm/day)','FontSize',9);
title('AM-65.5W-7.5S','fontsize',14);
text(0.5,11.4,'a1','FontSize',7,'FontWeight','bold');
grid on;
subplot1(2);
plot(xx,rain_africa,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('AF-23.5E-2.5S','fontsize',14);
text(0.5,7.5,'a2','FontSize',7,'FontWeight','bold');
grid on;
subplot1(3);
plot(xx,rain_asia,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('MA-113.5E-1.5N','fontsize',14);
text(0.5,9.5,'a3','FontSize',7,'FontWeight','bold');
grid on;
subplot1(4);
plot(xx,rain_america,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('US-83.5W-34.5N','fontsize',14);
text(0.5,14,'a4','FontSize',7,'FontWeight','bold');
grid on;
subplot1(5);
plot(xx,rain_europe,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'','','','','','','','','','','',''});
title('EU-4.5E-45.5N','fontsize',14);
text(0.5,5.5,'a5','FontSize',7,'FontWeight','bold');
grid on;
subplot1(6);
plot(xx,elai_amazon,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:12);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
ylabel('Leaf area index','FontSize',9);
text(0.5,5.9,'b1','FontSize',7,'FontWeight','bold');
grid on;
subplot1(7);
plot(xx,elai_africa,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:12);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,5.9,'b2','FontSize',7,'FontWeight','bold');
grid on;
subplot1(8);
plot(xx,elai_asia,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,5.35,'b3','FontSize',7,'FontWeight','bold');
grid on;
subplot1(9);
plot(xx,elai_america,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,4.75,'b4','FontSize',7,'FontWeight','bold');
grid on;
subplot1(10);
plot(xx,elai_europe,'k-','LineWidth',1);hold off;
set(gca,'Xlim',[0 11]);
set(gca,'XTick',0:1:11);
set(gca,'XTickLabel',{'1','','3','','5','','7','','9','','11',''});
xlabel('Month','FontSize',9);
text(0.5,2.93,'b5','FontSize',7,'FontWeight','bold');
grid on;
set(gcf,'PaperPositionMode','auto')
savefigure('Revised_Figure3_Seasonal_rain_lai.tif','PrintOptions','-dtiff', 'Size',[6.8 3],'LineWidth',1)
disp('done'); |
|