pressure_altitude (matlab)

最新推荐文章于 2022-12-12 15:17:01 发布

NAN_WHY

最新推荐文章于 2022-12-12 15:17:01 发布

阅读量513

点赞数

分类专栏： matlab

本文链接：https://blog.youkuaiyun.com/NAN_WHY/article/details/82194917

版权

matlab 专栏收录该内容

0 篇文章

订阅专栏

该博客主要围绕CMIP5数据进行气压高度计算。首先加载不同时期的模型数据和观测数据，对气压（psl）和温度（tas）数据进行网格重采样。接着将模型数据插值到站点，最后进行边际偏差校正，计算不同时期的气压高度，并保存校正后的数据。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

clear all
close all

% calculation of pressure altitude with CMIP5 data psl(p0) tas(t)
% 05/07/2018 By Nan Zhang

file=dir('F:\CMIP5 from HP');
%choose the resolution
lon1=0:1:360;
lat1=-90:1:90;
[long,latg]=meshgrid(lon1,lat1);

for a=9:11 %model name
    %% load model data (1976-2005)
    %load psl and regrid
    eval(['file_psl=dir(''F:\CMIP5 from HP\',file(a).name,'\historical\Atmos\psl'');']);
    eval(['lon=ncread(''',file_psl(3).name,''',''lon'');']);
    eval(['lat=ncread(''',file_psl(3).name,''',''lat'');']);
          lon=roundn(lon,-2);
        lat=roundn(lat,-2);

    eval(['psl11=ncread(''',file_psl(3).name,''',''psl'');']);
    eval(['psl12=ncread(''',file_psl(4).name,''',''psl'');']);
    eval(['psl13=ncread(''',file_psl(5).name,''',''psl'');']);
    psl1=cat(3,psl11,psl12,psl13);
    psl1=maskregrid(psl1,lon,lat,lon1,lat1);

    %load tas and regrid
    eval(['file_tas=dir(''F:\CMIP5 from HP\',file(a).name,'\historical\Atmos\tas'');']);
    eval(['tas11=ncread(''',file_tas(3).name,''',''tas'');']);
    eval(['tas12=ncread(''',file_tas(4).name,''',''tas'');']);
    eval(['tas13=ncread(''',file_tas(5).name,''',''tas'');']);
    tas1=cat(3,tas11,tas12,tas13);
    tas1=maskregrid(tas1,lon,lat,lon1,lat1);


    %% load model data (2021-2050)(2071-2100)
    %load psl and regrid
    eval(['file_psl=dir(''F:\CMIP5 from HP\',file(a).name,'\rcp85\Atmos\psl'');']);
    eval(['psl21=ncread(''',file_psl(3).name,''',''psl'');']);
    eval(['psl22=ncread(''',file_psl(4).name,''',''psl'');']);
    eval(['psl23=ncread(''',file_psl(5).name,''',''psl'');']);
    eval(['psl31=ncread(''',file_psl(6).name,''',''psl'');']);
    eval(['psl32=ncread(''',file_psl(7).name,''',''psl'');']);
    eval(['psl33=ncread(''',file_psl(8).name,''',''psl'');']);
    psl2=cat(3,psl21,psl22,psl23);
    psl3=cat(3,psl31,psl32,psl33);
    psl2=maskregrid(psl2,lon,lat,lon1,lat1);
    psl3=maskregrid(psl3,lon,lat,lon1,lat1);

    %load tas and regrid
    eval(['file_tas=dir(''F:\CMIP5 from HP\',file(a).name,'\rcp85\Atmos\tas'');']);
    eval(['tas21=ncread(''',file_tas(3).name,''',''tas'');']);
    eval(['tas22=ncread(''',file_tas(4).name,''',''tas'');']);
    eval(['tas23=ncread(''',file_tas(5).name,''',''tas'');']);
    eval(['tas31=ncread(''',file_tas(6).name,''',''tas'');']);
    eval(['tas32=ncread(''',file_tas(7).name,''',''tas'');']);
    eval(['tas33=ncread(''',file_tas(8).name,''',''tas'');']);
    tas2=cat(3,tas21,tas22,tas23);
    tas3=cat(3,tas31,tas32,tas33);
    tas2=maskregrid(tas2,lon,lat,lon1,lat1);
    tas3=maskregrid(tas3,lon,lat,lon1,lat1);


       %% load observational data
    station=xlsread('机场站点信息国外.xlsx');
    stationlon=station(:,1); %station longitude and latitude
    stationlat=station(:,3);
    h0=station(:,8); %airport altitude above sea level

    %% interpolate model data into stations(nearest neighborhood interpolation)
    TAS1=[];TAS2=[];TAS3=[];PSL1=[];PSL2=[];PSL3=[];
    %1976-2005
    for i=1:size(psl1,3) % the number of days during the 30 years
        tas1s=interp2(long,latg,tas1(:,:,i)',stationlon,stationlat,'nearest'); % bilinear interpolation day by day
        TAS1=[TAS1,tas1s]; % station*date
        psl1s=interp2(long,latg,psl1(:,:,i)',stationlon,stationlat,'nearest');
        PSL1=[PSL1,psl1s];

    end
    %2021-2050
    for i=1:size(tas2,3)
        tas2s=interp2(long,latg,tas2(:,:,i)',stationlon,stationlat,'nearest');
        TAS2=[TAS2,tas2s];
        psl2s=interp2(long,latg,psl2(:,:,i)',stationlon,stationlat,'nearest');
        PSL2=[PSL2,psl2s];

    end
    %2071-2100
    for i=1:size(tas3,3)
        tas3s=interp2(long,latg,tas3(:,:,i)',stationlon,stationlat,'nearest');
        TAS3=[TAS3,tas3s];
        psl3s=interp2(long,latg,psl3(:,:,i)',stationlon,stationlat,'nearest');
        PSL3=[PSL3,psl3s];

    end

    TAS1=TAS1';%date*station
    TAS2=TAS2';
    TAS3=TAS3';
    PSL1=PSL1';
    PSL2=PSL2';
    PSL3=PSL3';


    %% marginal bias correction (density altitude)
    % t(TAS)=temperature, K
    % p(P)=absolte pressure, Pa
    % psl(PSL)=sea level pressure, Pa
    % h0=height above sea level, m
    % hp(HP)=pressure altitude, m

    load('ObservationAirportsInternational.mat');
    sta1=STATION10(:,2);
    idd=find(mod(STATION10(:,2),10000)==229);
    sta2=sta1;
    sta2(idd)=[];

   P1=[];P2=[];P3=[];HP1C=[];HP2C=[];HP3C=[];TAS1C=[];TAS2C=[];TAS3C=[];% HP1C:纠偏以后的pressure altitude
    for j=1:18
        eval(['data=STATION',num2str(j),';']);
      if size(tas3,3)==10950 %judge the model whether contains the 29th day of leap month
          id=find(mod(data(:,2),10000)==229);
          data(id,:)=[];
          [c,ia]=setdiff(sta2,data(:,2));
      else [c,ia]=setdiff(sta1,data(:,2));
      end

         t=data(:,3); psl=data(:,4);
         site=find(t==9999.9 | psl==9999.9);
         t(site)=[];
         t=(t-32)/1.8+273.15; % convert from ℉ to K
         psl(site)=[];
         psl=psl*100; % convert from mb to Pa

         %ho,to calculate observational pressure altitude 1976-2005
        p=psl.*(t./(t+0.0065*h0(j))).^5.257;
        hp=(1-(p./101325).^0.190284)*145366.45*0.3048;

        %hs, to calculate simulated density altitude 1976-2005 for calibration
        PSL=PSL1;PSL(ia,:)=[];PSL(site,:)=[];PSL=PSL(:,j);
        TAS=TAS1;TAS(ia,:)=[];TAS(site,:)=[];TAS=TAS(:,j);
        P=PSL.*(TAS./(TAS+0.0065*h0(j))).^5.257;
        HP=(1-(P./101325).^0.190284)*145366.45*0.3048;

        %fs(1976-2005)
        TAS1s=TAS1(:,j);
        PSL1s=PSL1(:,j);
        P1s=PSL1s.*(TAS1s./(TAS1s+0.0065*h0(j))).^5.257;
        HP1s=(1-(P1s./101325).^0.190284)*145366.45*0.3048;
        P1=[P1,P1s];

        %fs(2021-2050)
        TAS2s=TAS2(:,j);
        PSL2s=PSL2(:,j);
        P2s=PSL2s.*(TAS2s./(TAS2s+0.0065*h0(j))).^5.257;
        HP2s=(1-(P2s./101325).^0.190284)*145366.45*0.3048;
        P2=[P2,P2s];

        %fs(2071-2100)
        TAS3s=TAS3(:,j);
        PSL3s=PSL3(:,j);
        P3s=PSL3s.*(TAS3s./(TAS3s+0.0065*h0(j))).^5.257;
        HP3s=(1-(P3s./101325).^0.190284)*145366.45*0.3048;
        P3=[P3,P3s];


        %mbc1976-2005
        HPFC1=mbc('temp',HP,hp,HP1s);
        HP1C=[HP1C,HPFC1];
        TASFC1=mbc('temp',TAS,t,TAS1s);
        TAS1C=[TAS1C,TASFC1];

        %mbc2021-2050
        HPFC2=mbc('temp',HP,hp,HP2s);
        HP2C=[HP2C,HPFC2];
        TASFC2=mbc('temp',TAS,t,TAS2s);
        TAS2C=[TAS2C,TASFC2];

        %mbc2071-2100
        HPFC3=mbc('temp',HP,hp,HP3s);
        HP3C=[HP3C,HPFC3];
        TASFC3=mbc('temp',TAS,t,TAS3s);
        TAS3C=[TAS3C,TASFC3];
    end
    eval(['save ',file(a).name,'.mat HP1C HP2C HP3C TAS1 TAS2 TAS3 TAS1C TAS2C TAS3C PSL1 PSL2 PSL3 P1 P2 P3 ']);
    % h0: 18 airports' height above sea level
end