Sample analysis code

Author: gaelforget

sample_analysis/example_IO.m

@@ -0,0 +1,53 @@
+function []=example_IO();
+%
+% EXAMPLE_IO computes and displays a standard deviation field
+%
+% stand-alone call: addpath gcmfaces/sample_analysis/; example_IO;
+%
+% needed input files:
+%   mkdir release1
+%   wget --recursive ftp://mit.ecco-group.org/ecco_for_las/version_4/release1/nctiles_climatology/ETAN
+%   mv mit.ecco-group.org/ecco_for_las/version_4/release1/nctiles_climatology release1/.
+
+gcmfaces_global;
+
+input_list_check('example_IO',nargin);
+
+%expected location:
+myenv.nctilesdir=fullfile(pwd,'/release1/nctiles_climatology/ETAN/');
+%if ETAN is not found then try old location:
+if ~isdir(myenv.nctilesdir);
+    %if not found then try old location:
+    tmpdir=fullfile(pwd,'/gcmfaces/sample_input/nctiles_climatology/ETAN/');
+    if isdir(tmpdir); myenv.nctilesdir=tmpdir; end;
+end;
+%if ETAN is still not found then issue warning and skip example_IO
+if ~isdir(myenv.nctilesdir);
+    help example_IO;
+    warning('example_IO requires release1/nctiles_climatology/ETAN that was not found ---> abort!');
+    return;
+end;
+
+%%%%%%%%%%%%%%%%%
+%load grid:
+%%%%%%%%%%%%%%%%%
+
+if isempty(mygrid);
+   grid_load;
+end;
+
+%%%%%%%%%%%
+%get field:
+%%%%%%%%%%%
+fld=read_nctiles([myenv.nctilesdir 'ETAN'],'ETAN');
+fld=std(fld,[],3); fld(find(fld==0))=NaN;
+cc=[0:0.1:1]*0.10;
+
+%%%%%%%%%%%%
+%plot field:
+%%%%%%%%%%%%
+if ~myenv.lessplot;
+    figureL; gcmfaces_sphere(fld,cc,[],[],1);
+end;
+
+

sample_analysis/example_budget.m

@@ -0,0 +1,133 @@
+function []=example_budget();
+%EXAMPLE_TRANSPORTS illustrates ECCO v4 budgets
+%
+%stand-alone call: addpath gcmfaces/sample_analysis/; example_budget;
+%
+%needed input files:
+%  wget --recursive ftp://mit.ecco-group.org/gforget/nctiles_budget_2d
+%  mv mit.ecco-group.org/gforget/nctiles_budget_2d sample_input/.
+%  rm -rf mit.ecco-group.org
+%
+%note: in the nctiles_budget_2d files from 03-Feb-2015 the sign convention 
+%   for trWtop and trWbot was positive downward; this sign convention was
+%   reversed in later versions for consistency with MITgcm convention.
+
+gcmfaces_global;
+
+if myenv.verbose>0;
+    gcmfaces_msg('===============================================');
+    gcmfaces_msg(['*** entering example_budget ' ...
+        'load budget terms from nctiles files ' ...
+        'and compute hemispheric budgets'],'');
+end;
+
+%%%%%%%%%%%%%%%%%
+%load grid:
+%%%%%%%%%%%%%%%%%
+
+%expected location:
+myenv.nctilesdir=fullfile('sample_input',filesep,'nctiles_budget_2d',filesep);
+%if nctiles_budget_2d is not found then try old location:
+if ~isdir(myenv.nctilesdir);
+    %if not found then try old location:
+    tmpdir=fullfile(myenv.gcmfaces_dir,'/sample_input/nctiles_budget_2d/');
+    if isdir(tmpdir); myenv.nctilesdir=tmpdir; end;
+end;
+%if nctiles_budget_2d is still not found then issue warning and skip example_budget
+if ~isdir(myenv.nctilesdir);
+    diags=[];
+    help example_budget;
+    warning(['skipping example_budget (missing ' myenv.nctilesdir ')']);
+    return;
+end;
+
+dirName=fullfile(myenv.nctilesdir,filesep,'budgMo',filesep);
+
+if ~isdir(dirName);
+    help gcmfaces_demo;
+    warning(['skipping example_budget (missing ' dirName ')']);
+    return;
+end;
+
+if isempty(which('ncload'));
+    warning(['skipping example_budget (missing ncload that is part of MITprof)']);
+    return;
+end;
+
+%%%%%%%%%%%%%%%%%%
+%main computation:
+%%%%%%%%%%%%%%%%%%
+
+%load grid:
+if isempty(mygrid);
+   grid_load;
+end;
+
+%select budget of interest
+nameBudg='budgMo';
+
+%load budget terms
+fileName=[myenv.nctilesdir nameBudg filesep nameBudg(1:6)];
+tend=read_nctiles(fileName,'tend');
+trU=read_nctiles(fileName,'trU');
+trV=read_nctiles(fileName,'trV');
+trWtop=read_nctiles(fileName,'trWtop');
+trWbot=read_nctiles(fileName,'trWbot');
+
+%load dt (time increments) vector (not a gcmfaces object)
+ncload([fileName '.0001.nc'],'dt');
+
+%get budget descitption and units
+ncid = netcdf.open([fileName '.0001.nc'],'NC_NOWRITE');
+descr = netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description');
+varid = netcdf.inqVarID(ncid,'tend');
+units = netcdf.getAtt(ncid,varid,'units');
+netcdf.close(ncid);
+
+if myenv.verbose>0; gcmfaces_msg(['* display Northern Hemisphere budget for ' nameBudg ' (' descr ')']);end;
+
+%define northern hemisphere as domain of integration
+nameMask='Northern Hemisphere';
+mask=mygrid.mskC(:,:,1).*(mygrid.YC>0); 
+areaMask=mygrid.RAC.*mask;
+
+%edit plot title accordingly
+tmp1=strfind(descr,'-- ECCO v4');
+descr=[descr(1:tmp1-1) 'for: ' nameMask];
+
+%compute northern hemisphere integrals
+budg.tend=NaN*dt;
+budg.hconv=NaN*dt;
+budg.zconv=NaN*dt;
+for tt=1:length(dt);
+    %compute flux convergence
+    hconv=calc_UV_conv(trU(:,:,tt),trV(:,:,tt));
+    zconv=trWtop(:,:,tt)-trWbot(:,:,tt);
+    %compute sum over domain
+    budg.tend(tt)=nansum(tend(:,:,tt).*mask)/nansum(areaMask);
+    budg.hconv(tt)=nansum(hconv.*mask)/nansum(areaMask);
+    budg.zconv(tt)=nansum(zconv.*mask)/nansum(areaMask);
+end;
+
+%display result
+figureL;
+subplot(3,1,1); set(gca,'FontSize',12);
+plot(cumsum(dt.*budg.tend));
+grid on; xlabel('month'); ylabel([units ' x s']);
+legend('content anomaly');
+subplot(3,1,2); set(gca,'FontSize',12);
+plot(cumsum(dt.*budg.tend)); hold on;
+plot(cumsum(dt.*budg.hconv),'r'); 
+plot(cumsum(dt.*budg.zconv),'g');
+grid on; xlabel('month'); ylabel([units ' x s']);
+legend('content anomaly','horizontal convergence','vertical convergence');
+subplot(3,1,3); set(gca,'FontSize',12);
+plot(cumsum(dt.*(budg.tend-budg.hconv-budg.zconv))); 
+grid on; xlabel('month'); ylabel([units ' x s']);
+legend('budget residual');
+
+if myenv.verbose>0;
+    gcmfaces_msg('*** leaving example_budget ','');
+    gcmfaces_msg('===============================================');
+end;
+

sample_analysis/example_display.m

@@ -0,0 +1,128 @@
+function []=example_display(method);
+% EXAMPLE_DISPLAY(method) illustrate various options to display 2D fields.
+%
+%  method (optional) selects the plotting method
+%     1: plot a gcmfaces object face by face (as in Fig. C1) of doi:10.5194/gmd-8-3071-2015).
+%     2: use convert2array via qwckplot (a quick but crude display option).
+%     3: use convert2pcol then pcolor (for a display in lat-lon coordinates).
+%     4: use the m_map_gcmfaces front end to m_map (https://www.eoas.ubc.ca/~rich/map.html).
+%     5: use gcmfaces_sphere to display fields on a sphere.
+%     6: use stretched vertical coordinate.
+%  If method is not specified then it is set to [3:5] by default.
+%
+% Example: addpath gcmfaces/sample_analysis/; example_display;
+
+gcmfaces_global;
+
+input_list_check('example_display',nargin);
+
+if isempty(whos('method')); method=[3:5]; end;
+
+if myenv.verbose>0;
+    gcmfaces_msg('===============================================');
+    gcmfaces_msg(['*** entering example_display: displays a gridded ' ...
+        'field of ocean depth (gcmfaces object) in geographic coordinates ' ...
+        'using pcolor, in various projections (if m_map is in Matlab path), ' ...
+        'and on a sphere.'],'');
+end;
+
+%%%%%%%%%%%%%%%%%
+%load grid:
+%%%%%%%%%%%%%%%%%
+
+if isempty(mygrid);
+   grid_load; 
+end;
+nF=mygrid.nFaces;
+
+%%%%%%%%%%%
+%get field:
+%%%%%%%%%%%
+
+fld=mygrid.Depth; fld(fld==0)=NaN;
+cc=[[0:0.05:0.5] [0.6 0.75 1 1.25]]*1e4; myCmap='gray';
+
+%%%%%%%%%%%%
+%plot field:
+%%%%%%%%%%%%
+
+if sum(ismember(method,1));
+    if myenv.verbose>0; gcmfaces_msg('* gcmfaces format display -- face by face.'); end;
+    if nF==1;
+        figureL; imagescnan(fld{1}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+    elseif nF==5;
+        figureL;
+        subplot(3,3,7); imagescnan(fld{1}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+        tmp1=axis; tmp2=text(tmp1(2)/2,tmp1(4)/2,'1','FontSize',32,'Color','r','Rotation',0);
+        subplot(3,3,8); imagescnan(fld{2}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+        tmp1=axis; tmp2=text(tmp1(2)/2,tmp1(4)/2,'2','FontSize',32,'Color','r','Rotation',0);
+        subplot(3,3,5); imagescnan(fld{3}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+        tmp1=axis; tmp2=text(tmp1(2)/2,tmp1(4)/2,'3','FontSize',32,'Color','r','Rotation',0);
+        subplot(3,3,6); imagescnan(fld{4}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+        tmp1=axis; tmp2=text(tmp1(2)/2,tmp1(4)/2,'4','FontSize',32,'Color','r','Rotation',0);
+        subplot(3,3,3); imagescnan(fld{5}','nancolor',[1 1 1]*0.8); axis xy; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb);
+        tmp1=axis; tmp2=text(tmp1(2)/2,tmp1(4)/2,'5','FontSize',32,'Color','r','Rotation',0);
+    title('display face by face');
+    elseif nF==6;
+        error('face by face plot for cude : not yet implemented');
+    end;
+end;
+
+if sum(ismember(method,2));
+    if myenv.verbose>0; gcmfaces_msg('* array format display -- all faces concatenated.'); end;
+    figureL; qwckplot(fld); cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); %delete(cb);
+    title('display using qwckplot');
+    for ff=1:mygrid.nFaces;
+        tmp0=0*mygrid.XC;
+        tmp1=round(size(tmp0{ff})/2);
+        tmp0{ff}(tmp1(1),tmp1(2))=1;
+        tmp0=convert2array(tmp0); [ii,jj]=find(tmp0==1);
+        if ff==3; ang=90; elseif ff>3; ang=-90; else; ang=0; end;
+        hold on; text(ii,jj,num2str(ff),'FontSize',32,'Color','r','Rotation',ang);
+    end;
+end;
+
+if sum(ismember(method,3));
+    if myenv.verbose>0; gcmfaces_msg('* geographical display -- using pcolor directly.'); end;
+    figureL;
+    [X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,fld); pcolor(X,Y,FLD);
+    if ~isempty(find(X>359)); axis([0 360 -90 90]); else; axis([-180 180 -90 90]); end;
+    shading flat; cb=gcmfaces_cmap_cbar(cc,{'myCmap',myCmap}); delete(cb); 
+    xlabel('longitude'); ylabel('latitude'); 
+    title('display using convert2pcol');
+end;
+
+if sum(ismember(method,4));
+    if myenv.verbose>0; gcmfaces_msg('* geographical display -- using m_map_gcmfaces.'); end;
+    if ~isempty(which('m_proj'));
+        figureL; m_map_gcmfaces(fld,0,{'myCaxis',cc},{'myCmap',myCmap});
+        aa=get(gcf,'Children'); axes(aa(4));
+        title('display using m_map_gcmfaces','Interpreter','none');
+    elseif myenv.verbose;
+        fprintf('  > To use m_map_gcmfaces, please add m_map to your Matlab path\n');
+    end;
+end;
+
+if sum(ismember(method,5));
+    if myenv.verbose>0; gcmfaces_msg('* geographical display -- on a sphere.'); end;
+    figureL; gcmfaces_sphere(fld,cc,[],'N',3);
+    title('display using gcmfaces_sphere','Interpreter','none');
+end;
+
+%test case for depthStretchPlot:
+if sum(ismember(method,6));
+    if myenv.verbose>0; gcmfaces_msg('* section display -- using strecthed vertical coord.'); end;
+    x=ones(length(mygrid.RC),1)*[1:200]; z=mygrid.RC*ones(1,200); c=sin(z/2000*pi).*cos(x/50*pi);
+    figureL;
+    subplot(1,2,1); set(gca,'FontSize',16);
+    pcolor(x,z,c); shading flat; title('standard depth display');
+    subplot(1,2,2); set(gca,'FontSize',16);
+    depthStretchPlot('pcolor',{x,z,c}); shading flat; 
+    title('stretched depth display');
+end;
+
+if myenv.verbose>0;
+    gcmfaces_msg('*** leaving example_display');
+    gcmfaces_msg('===============================================');
+end;
+