ECCO v4r5处理

起风了my

于 2024-07-24 17:44:19 发布

阅读量296

点赞数 20

文章标签： matlab

本文链接：https://blog.csdn.net/weixin_42510703/article/details/140649040

版权

参考资料：

文章：https://egusphere.copernicus.org/preprints/2024/egusphere-2024-727/egusphere-2024-727.pdf
ECCOv4r4数据用户指南：https://ecco-group.org/user-guide-v4r4.htm
ECCOv4r5 资料pdf1：https://ecco-group.org/docs/ecco_annual_mtg23_day2_01_fukumori.pdf
ECCOv4r5出图pdf2：https://ecco.jpl.nasa.gov/drive/files/Version4/Release5/doc/v4r5_overview_plots.pdf
ECCO模拟与观测对比结果展示 Model/Observation Comparison in Recent ECCOs: Sea ice Fields in Southern Ocean:https://ecco.odyseallc.net/docs/03_01_Zhang_ecco_Seaice_hong_zhang.pdf
ECCO support（不会的可以发邮件请教）:https://mailman.mit.edu/mailman/listinfo/ecco-support
gcmfaces库下载地址： https://gcmfaces.readthedocs.io/en/latest/prep_install.html
数据下载网址：https://engaging-web.mit.edu/~gforget/harbor/version_4/release2/
MITgcm下载网址：http://wwwcvs.mitgcm.org/viewvc/MITgcm/MITgcm/utils/matlab/

本文目的：看一下ECCOv4r5的海冰面积模拟情况。
主要解决的问题：如何用matlab打开.data数据，并且成功出图看Sea ice area的情况。

第一步：下载数据和代码：

参考资料1文章对应出图代码下载地址：https://zenodo.org/records/10958934
本文具体参照的代码是ECCOv4r5/Seaice_eccov4.m
在这里插入图片描述
ECCOv4r5数据下载网址：https://zenodo.org/records/10930853
下载我们需要的SIarea_mon_mean 海冰面积月数据，为了实现作者代码，还需要下载SIhef冰架数据

MITgcm网址中需下载的cs_grid文件夹：

第二步：
根据代码需要的文件，整理好文件夹的结构，并下载需要的库，例如m_map，cmocean等。
代码所需的文件
然后将其中的文件夹改成自己电脑的路径
在这里插入图片描述
我文件夹的结构如下图所示：

第三步：用matlab读取ECCOv4r5数据并出图

% Let`s make figures of cross section about temperature/salinity/potential density every 10 days for 3 months.
% cd /Volumes/Samsung_T5/ECCO_v4r5/Version4/Release5/diags_monthly
clc;close;clear all;
plotA=1;
plotH=0;

%path
addpath C:/STUDY/MITgcm-master/utils/matlab/cs_grid/read_cs/
addpath C:/STUDY/MITgcm-master/utils/matlab/cs_grid/
addpath C:/STUDY/MITgcm-master/utils/matlab/
addpath C:/STUDY/MITgcm-master/utils/matlab/cs_grid/latloncap/
% addpath C:/STUDY/MITgcm_master/matlab/m_map/
addpath C:/STUDY/cmocean/

fnamx=['C:/STUDY/grid/XG.data'];
fnamy=['C:/STUDY/grid/YC.data'];
fnamdyg=['C:/STUDY/grid/RAC.data']; %Caution RAC
fnamhFacW=['C:/STUDY/grid/hFacC.data'];
RC=rdmds('C:/STUDY/grid/RC');
DRF=rdmds('C:/STUDY/grid/DRF');
%RAC=rdmds('./grid/RAC');
fc=[1 2 3 4 5];
%这个时间表示每一年的每个月；
time=[732,   1428,  2172,  2892,  3636,  4356,  5100,  5844,  6564,  7308,  8028,  8772, ...
      9516, 10188, 10932, 11652, 12396, 13116, 13860, 14604, 15324, 16068, 16788, 17532, ...
     18276, 18948, 19692, 20412, 21156, 21876, 22620, 23364, 24084, 24828, 25548, 26292, ...
     27036, 27708, 28452, 29172, 29916, 30636, 31380, 32124, 32844, 33588, 34308, 35052, ...
     35796, 36492, 37236, 37956, 38700, 39420, 40164, 40908, 41628, 42372, 43092, 43836, ...
     44580, 45252, 45996, 46716, 47460, 48180, 48924, 49668, 50388, 51132, 51852, 52596, ...
     53340, 54012, 54756, 55476, 56220, 56940, 57684, 58428, 59148, 59892, 60612, 61356, ...
     62100, 62772, 63516, 64236, 64980, 65700, 66444, 67188, 67908, 68652, 69372, 70116, ...
     70860, 71556, 72300, 73020, 73764, 74484, 75228, 75972, 76692, 77436, 78156, 78900, ...
     79644, 80316, 81060, 81780, 82524, 83244, 83988, 84732, 85452, 86196, 86916, 87660, ...
     88404, 89076, 89820, 90540, 91284, 92004, 92748, 93492, 94212, 94956, 95676, 96420, ...
     97164, 97836, 98580, 99300,100044,100764,101508,102252,102972,103716,104436,105180, ...
    105924,106620,107364,108084,108828,109548,110292,111036,111756,112500,113220,113964, ...
    114708,115380,116124,116844,117588,118308,119052,119796,120516,121260,121980,122724, ...
    123468,124140,124884,125604,126348,127068,127812,128556,129276,130020,130740,131484, ...
    132228,132900,133644,134364,135108,135828,136572,137316,138036,138780,139500,140244, ...
    140988,141684,142428,143148,143892,144612,145356,146100,146820,147564,148284,149028, ...
    149772,150444,151188,151908,152652,153372,154116,154860,155580,156324,157044,157788, ...
    158532,159204,159948,160668,161412,162132,162876,163620,164340,165084,165804,166548, ...
    167292,167964,168708,169428,170172,170892,171636,172380,173100,173844,174564,175308, ...
    176052,176748,177492,178212,178956,179676,180420,181164,181884,182628,183348,184092, ...
    184836,185508,186252,186972,187716,188436,189180,189924,190644,191388,192108,192852, ...
    193596,194268,195012,195732,196476,197196,197940,198684,199404,200148,200868,201612, ...
    202356,203028,203772,204492,205236,205956,206700,207444,208164,208908,209628,210372, ...
    211116,211812,212556,213276,214020,214740,215484,216228,216948,217692,218412,219156, ...
    219900,220572,221316,222036,222780,223500,224244,224988,225708,226452,227172,227916, ...
    228660,229332,230076,230796,231540,232260,233004,233748,234468,235212,235932,236676, ...
    237420,238092,238836,239556,240300,241020,241764,242508,243228,243972,244692,245424];

tt=1; %237420
nx=90;

 f=1; fldx1(:,:)=read_llc_fkij(fnamx,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=2; fldx2(:,:)=read_llc_fkij(fnamx,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=3; fldx3(:,:)=read_llc_fkij(fnamx,nx,fc(f),1,1:nx,1:  nx,'real*4');
 f=4; fldx4(:,:)=read_llc_fkij(fnamx,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=5; fldx5(:,:)=read_llc_fkij(fnamx,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 
 f=1; fldy1(:,:)=read_llc_fkij(fnamy,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=2; fldy2(:,:)=read_llc_fkij(fnamy,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=3; fldy3(:,:)=read_llc_fkij(fnamy,nx,fc(f),1,1:nx,1:  nx,'real*4');
 f=4; fldy4(:,:)=read_llc_fkij(fnamy,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=5; fldy5(:,:)=read_llc_fkij(fnamy,nx,fc(f),1,1:nx,1:3*nx,'real*4');


 f=1; flddyg1(:,:)=read_llc_fkij(fnamdyg,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=2; flddyg2(:,:)=read_llc_fkij(fnamdyg,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=3; flddyg3(:,:)=read_llc_fkij(fnamdyg,nx,fc(f),1,1:nx,1:  nx,'real*4');
 f=4; flddyg4(:,:)=read_llc_fkij(fnamdyg,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 f=5; flddyg5(:,:)=read_llc_fkij(fnamdyg,nx,fc(f),1,1:nx,1:3*nx,'real*4');
 
 for k=1:50 
 f=1; fldhFacW1(:,:,k)=read_llc_fkij(fnamhFacW,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=2; fldhFacW2(:,:,k)=read_llc_fkij(fnamhFacW,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=3; fldhFacW3(:,:,k)=read_llc_fkij(fnamhFacW,nx,fc(f),k,1:nx,1:  nx,'real*4');
 f=4; fldhFacW4(:,:,k)=read_llc_fkij(fnamhFacW,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=5; fldhFacW5(:,:,k)=read_llc_fkij(fnamhFacW,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 end

XC(1:90,:)=fldx1; XC(91:180,:)=fldx2; XC(181:270,:)=fldx4; XC(271:360,:)=fldx5;  
YC(1:90,:)=fldy1; YC(91:180,:)=fldy2; YC(181:270,:)=fldy4; YC(271:360,:)=fldy5;  
RAC(1:90,:)=flddyg1; RAC(91:180,:)=flddyg2; RAC(181:270,:)=flddyg4; RAC(271:360,:)=flddyg5; 
hFacW(1:90,:,:)=fldhFacW1; hFacW(91:180,:,:)=fldhFacW2; hFacW(181:270,:,:)=fldhFacW4;  hFacW(271:360,:,:)=fldhFacW5; 

%tset=[12, 60, 120, 336]

ccc=0;
%1992-2019年数据，共28年
%tttt表示一年12个月
for yr=1:28
for tttt=1:12
    ccc=ccc+1
tt=time(ccc)
if(ccc<2); 
  fnamU=['C:/STUDY/SIarea_mon_mean/SIarea_mon_mean.0000000' num2str(time(ccc)) '.data'];
elseif(ccc<14);
  fnamU=['C:/STUDY/SIarea_mon_mean/SIarea_mon_mean.000000' num2str(time(ccc)) '.data'];
elseif(ccc<137);  
  fnamU=['C:/STUDY/SIarea_mon_mean/SIarea_mon_mean.00000' num2str(time(ccc)) '.data'];
else
  fnamU=['C:/STUDY/SIarea_mon_mean/SIarea_mon_mean.0000' num2str(time(ccc)) '.data'];
end

if(ccc<2); 
  fnamV=['C:/STUDY/SIheff_mon_mean/SIheff_mon_mean.0000000' num2str(time(ccc)) '.data'];
elseif(ccc<14);
  fnamV=['C:/STUDY/SIheff_mon_mean/SIheff_mon_mean.000000' num2str(time(ccc)) '.data'];
elseif(ccc<137);  
  fnamV=['C:/STUDY/SIheff_mon_mean/SIheff_mon_mean.00000' num2str(time(ccc)) '.data'];
else
  fnamV=['C:/STUDY/SIheff_mon_mean/SIheff_mon_mean.0000' num2str(time(ccc)) '.data'];
end
 
for k=1:1
 f=1; fldT1(:,:,k)=read_llc_fkij(fnamU,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=2; fldT2(:,:,k)=read_llc_fkij(fnamU,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=3; fldT3(:,:,k)=read_llc_fkij(fnamU,nx,fc(f),k,1:nx,1:  nx,'real*4');
 f=4; fldT4(:,:,k)=read_llc_fkij(fnamU,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=5; fldT5(:,:,k)=read_llc_fkij(fnamU,nx,fc(f),k,1:nx,1:3*nx,'real*4');     
     
 f=1; fldS1(:,:,k)=read_llc_fkij(fnamV,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=2; fldS2(:,:,k)=read_llc_fkij(fnamV,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=3; fldS3(:,:,k)=read_llc_fkij(fnamV,nx,fc(f),k,1:nx,1:  nx,'real*4');
 f=4; fldS4(:,:,k)=read_llc_fkij(fnamV,nx,fc(f),k,1:nx,1:3*nx,'real*4');
 f=5; fldS5(:,:,k)=read_llc_fkij(fnamV,nx,fc(f),k,1:nx,1:3*nx,'real*4'); 
end
%T表示海冰面积，S表示冰架
T(1:90,:,yr,tttt)=fldT1; T(91:180,:,yr,tttt)=fldT2; T(181:270,:,yr,tttt)=fldT4; T(271:360,:,yr,tttt)=fldT5;
S(1:90,:,yr,tttt)=fldS1; S(91:180,:,yr,tttt)=fldS2; S(181:270,:,yr,tttt)=fldS4; S(271:360,:,yr,tttt)=fldS5;

end
end

AA2=T; 
HH2=S;

%mypcolor(streamfunction); caxis([-40,40]);


for i=1:360
 for j=1:270
  for k=1:1
    if(hFacW(i,j,k)==0); AA2(i,j,:,:)=NaN; HH2(i,j,:,:)=NaN; end
  end
 end
end


for i=1:360
 for j=1:270
MARCH_AREA(i,j)=mean(AA2(i,j,:,3)); %1992-2018
MARCH_HEFF(i,j)=mean(HH2(i,j,:,3)); %1992-2018
SEP_AREA(i,j)=mean(AA2(i,j,:,9)); %1992-2018
SEP_HEFF(i,j)=mean(HH2(i,j,:,9)); %1992-2018
 end
end

% if(plotA)
% subplot(1,2,1)
% SF_t=squeeze(MARCH_AREA(:,:)); 
% %convert2spherical_cordinate
% latn=102;
% 
% XC_t=XC(:,10:latn); XC_t(361,:)=XC(1,10:latn);
% YC_t=YC(:,10:latn); YC_t(361,:)=YC(1,10:latn);
% SF_t2=SF_t(:,10:latn); SF_t2(361,:)=SF_t(1,10:latn);
% XX2=cos(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% YY2=sin(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% h=pcolor(XX2,YY2,SF_t2(:,:)); 
% caxis([0., 1.]); colorbar; cmocean('ice');
% set(h, 'EdgeColor', 'none');
% set(gca,'fontsize',16)
% set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]);
% axis square
% 
% 
% subplot(1,2,2)
% SF_t=squeeze(SEP_AREA(:,:)); 
% %convert2spherical_cordinate
% latn=102;
% 
% XC_t=XC(:,10:latn); XC_t(361,:)=XC(1,10:latn);
% YC_t=YC(:,10:latn); YC_t(361,:)=YC(1,10:latn);
% SF_t2=SF_t(:,10:latn); SF_t2(361,:)=SF_t(1,10:latn);
% XX2=cos(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% YY2=sin(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% h=pcolor(XX2,YY2,SF_t2(:,:)); 
% caxis([0., 1.]); colorbar; cmocean('ice');
% set(h, 'EdgeColor', 'none');
% set(gca,'fontsize',16)
% set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]);
% axis square
% 
% end

% save ECCOv4r5_seaice.mat

%冰架
% {
% if(plots);
% for tts=1:4
% subplot(1,4,tts)
% ttt=tts;
% if(tts==1); SF_t=squeeze(S(:,:,24,tts)); end; %552 m
% if(tts>1); SF_t=squeeze(S(:,:,24,tts)-S(:,:,24,1)); end; %552 m
% % convert2spherical_cordinate
% latn=102;
% 
% XC_t=XC(:,10:latn); XC_t(361,:)=XC(1,10:latn);
% YC_t=YC(:,10:latn); YC_t(361,:)=YC(1,10:latn);
% SF_t2=SF_t(:,10:latn); SF_t2(361,:)=SF_t(1,10:latn);
% XX2=cos(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% YY2=sin(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% h=pcolor(XX2,YY2,SF_t2(:,:)); 
% if(tts==1); caxis([34.6 34.75]); colorbar; cmocean('haline'); end;
% if(tts>1); caxis([-0.1, 0.1]); colorbar; cmocean('balance'); end;
% set(h, 'EdgeColor', 'none');
% set(gca,'fontsize',16)
% set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]);
% axis square
% end
% end
% }



% %———————开始对1月份的海冰面积进行研究
JAN_AREA = AA2(:,:,:,1);
% %JAN_AREA的维度是360x270x28，其中28是年份数量，1992-2019年
% JAN_AREA
for k=1:28
subplot(4,7,k)
SF_t=squeeze(JAN_AREA(:,:,k)); 
%convert2spherical_cordinate
latn=102;

XC_t=XC(:,10:latn); XC_t(361,:)=XC(1,10:latn);
YC_t=YC(:,10:latn); YC_t(361,:)=YC(1,10:latn);
SF_t2=SF_t(:,10:latn); SF_t2(361,:)=SF_t(1,10:latn);
XX2=cos(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
YY2=sin(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
h=pcolor(XX2,YY2,SF_t2(:,:)); 
title(num2str(1991+k))
caxis([0., 1.]); colorbar; cmocean('ice');
set(h, 'EdgeColor', 'none');
set(gca,'fontsize',16)
set(gca,'YTick',0:0.2:1)
set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]);
axis square
end



% 
% subplot(1,2,2)
% SF_t=squeeze(SEP_AREA(:,:)); 
% %convert2spherical_cordinate
% latn=102;
% 
% XC_t=XC(:,10:latn); XC_t(361,:)=XC(1,10:latn);
% YC_t=YC(:,10:latn); YC_t(361,:)=YC(1,10:latn);
% SF_t2=SF_t(:,10:latn); SF_t2(361,:)=SF_t(1,10:latn);
% XX2=cos(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% YY2=sin(-XC_t*(2*pi/360)).*cos(YC_t*(2*pi/360));
% h=pcolor(XX2,YY2,SF_t2(:,:)); 
% caxis([0., 1.]); colorbar; cmocean('ice');
% set(h, 'EdgeColor', 'none');
% set(gca,'fontsize',16)
% set(gca,'YTickLabel',[]); set(gca,'XTickLabel',[]);
% axis square