| 1 |
function G = load_grid(varargin); |
| 2 |
% function G = load_grid(dirName,option,nyAxis,nxAxis); |
| 3 |
% |
| 4 |
% load the MITgcm output grid-files into one structure array |
| 5 |
% arguments: |
| 6 |
% dirName = directory where grid-files are |
| 7 |
% optional arguments: |
| 8 |
% ncdf = last digit of integer "option": |
| 9 |
% ncdf = 0,2 : read binary (MDSIO) files using rdmds (default: ncdf=0) |
| 10 |
% ncdf = 1,3 : read NetCDF (MNC) files using rdmnc |
| 11 |
% ncdf = 2,3 : as 0,1 + print elapsed-time for loading files |
| 12 |
% option >= 10 : only read in Horiz.Grid spacing (without hFac) |
| 13 |
% option >= 20 : only read in Verti.Grid spacing (without hFac) |
| 14 |
% nxAxis = number of spacing in 1.D x-axis xAxC (default: nxAxis = grid-Nx) |
| 15 |
|
| 16 |
% $Header: /u/gcmpack/MITgcm/utils/matlab/load_grid.m,v 1.1 2012/10/04 14:02:42 jmc Exp $ |
| 17 |
% $Name: $ |
| 18 |
|
| 19 |
if nargin == 0 |
| 20 |
rDir = '.'; |
| 21 |
else |
| 22 |
rDir = varargin{1}; |
| 23 |
end |
| 24 |
if nargin < 2, |
| 25 |
ktyp=0; |
| 26 |
ncdf=0; |
| 27 |
else |
| 28 |
ncdf= varargin{2}; |
| 29 |
ktyp=floor(ncdf/10); ncdf=rem(ncdf,10); |
| 30 |
end |
| 31 |
|
| 32 |
if ncdf > 1, TimeT0=clock; end |
| 33 |
|
| 34 |
if rem(ncdf,2) == 0, |
| 35 |
%- load MDSIO grid files : |
| 36 |
if ktyp < 2, |
| 37 |
xC=rdmds(fullfile(rDir,'XC')); |
| 38 |
yC=rdmds(fullfile(rDir,'YC')); |
| 39 |
xG=rdmds(fullfile(rDir,'XG')); |
| 40 |
yG=rdmds(fullfile(rDir,'YG')); |
| 41 |
|
| 42 |
dXc=rdmds(fullfile(rDir,'DXC')); |
| 43 |
dYc=rdmds(fullfile(rDir,'DYC')); |
| 44 |
dXg=rdmds(fullfile(rDir,'DXG')); |
| 45 |
dYg=rdmds(fullfile(rDir,'DYG')); |
| 46 |
|
| 47 |
rAc=rdmds(fullfile(rDir,'RAC')); |
| 48 |
rAw=rdmds(fullfile(rDir,'RAW')); |
| 49 |
rAs=rdmds(fullfile(rDir,'RAS')); |
| 50 |
rAz=rdmds(fullfile(rDir,'RAZ')); |
| 51 |
|
| 52 |
%- load grid orientation relative to West-East / South-North dir: |
| 53 |
if isempty(dir(fullfile(rDir,'AngleCS.*'))), |
| 54 |
fprintf(' no grid orientation (Cos & Sin) file to load '); |
| 55 |
csAngle= ones(size(rAc)); |
| 56 |
snAngle=zeros(size(rAc)); |
| 57 |
fprintf(' => set COS=1,SIN=0\n'); |
| 58 |
else |
| 59 |
csAngle=rdmds(fullfile(rDir,'AngleCS')); |
| 60 |
snAngle=rdmds(fullfile(rDir,'AngleSN')); |
| 61 |
end |
| 62 |
dims = size(rAc); |
| 63 |
end |
| 64 |
|
| 65 |
if rem(ktyp,2) == 0, |
| 66 |
rC=rdmds(fullfile(rDir,'RC')); rC=squeeze(rC); |
| 67 |
rF=rdmds(fullfile(rDir,'RF')); rF=squeeze(rF); |
| 68 |
dRc=rdmds(fullfile(rDir,'DRC')); dRc=squeeze(dRc); |
| 69 |
dRf=rdmds(fullfile(rDir,'DRF')); dRf=squeeze(dRf); |
| 70 |
dims = size(rC); |
| 71 |
end |
| 72 |
|
| 73 |
if ktyp == 0, |
| 74 |
%- define domain: |
| 75 |
hFacC=rdmds(fullfile(rDir,'hFacC')); |
| 76 |
hFacW=rdmds(fullfile(rDir,'hFacW')); |
| 77 |
hFacS=rdmds(fullfile(rDir,'hFacS')); |
| 78 |
depth=rdmds(fullfile(rDir,'Depth')); |
| 79 |
dims = size(hFacC); |
| 80 |
end |
| 81 |
|
| 82 |
if ncdf > 1, TimeT1=clock; end |
| 83 |
|
| 84 |
else |
| 85 |
%- load NetCDF grid files : |
| 86 |
%S=rdmnc([rDir,'grid.*.nc']); |
| 87 |
S=rdmnc(fullfile(rDir,'grid.*.nc')); |
| 88 |
if ncdf > 1, TimeT1=clock; end |
| 89 |
|
| 90 |
%-- |
| 91 |
xC=S.XC; |
| 92 |
yC=S.YC; |
| 93 |
xG=S.XG(1:end-1,1:end-1); |
| 94 |
yG=S.YG(1:end-1,1:end-1); |
| 95 |
rC=S.RC'; |
| 96 |
rF=S.RF'; |
| 97 |
dXc=S.dxC(1:end-1,:); |
| 98 |
dYc=S.dyC(:,1:end-1); |
| 99 |
dXg=S.dxG(:,1:end-1); |
| 100 |
dYg=S.dyG(1:end-1,:); |
| 101 |
dRc=S.drC'; |
| 102 |
dRf=S.drF'; |
| 103 |
rAc=S.rA; |
| 104 |
rAw=S.rAw(1:end-1,:); |
| 105 |
rAs=S.rAs(:,1:end-1); |
| 106 |
rAz=S.rAz(1:end-1,1:end-1); |
| 107 |
if isfield(S,'AngleCS') & isfield(S,'AngleSN'), |
| 108 |
csAngle=S.AngleCS; |
| 109 |
snAngle=S.AngleSN; |
| 110 |
else |
| 111 |
fprintf(' no grid orientation (Cos & Sin) in grid file '); |
| 112 |
csAngle= ones(size(rAc)); |
| 113 |
snAngle=zeros(size(rAc)); |
| 114 |
fprintf(' => set COS=1,SIN=0\n'); |
| 115 |
end |
| 116 |
hFacC=S.HFacC; |
| 117 |
hFacW=S.HFacW(1:end-1,:,:); |
| 118 |
hFacS=S.HFacS(:,1:end-1,:); |
| 119 |
depth=S.Depth; |
| 120 |
dims = size(hFacC); |
| 121 |
end |
| 122 |
|
| 123 |
if ncdf > 1, fprintf(' (took %6.4f s)\n',etime(TimeT1,TimeT0)); end |
| 124 |
|
| 125 |
if length(dims) == 1, dims(2)=1; end |
| 126 |
if length(dims) == 2, dims(3)=1; end |
| 127 |
|
| 128 |
%- 1-D axis: |
| 129 |
if nargin < 3, |
| 130 |
%- yAxis |
| 131 |
ny=dims(2); |
| 132 |
yAxC=yC(1,:)'; |
| 133 |
if rem(ncdf,2) == 0, |
| 134 |
yAxV=[yG(1,:) 2*yC(1,end)-yG(1,end)]'; |
| 135 |
else |
| 136 |
yAxV=S.YG(1,:)'; |
| 137 |
end |
| 138 |
else |
| 139 |
ny=varargin{3}; |
| 140 |
dyAx=(max(yG(:))-min(yG(:)))/ny; |
| 141 |
yAxV=min(yG(:))+[0:ny]'*dyAx; |
| 142 |
yAxC=(yAxV(1:ny)+yAxV(2:ny+1))/2; |
| 143 |
end |
| 144 |
if nargin < 4, |
| 145 |
%- xAxis |
| 146 |
nx=dims(1); |
| 147 |
xAxC=xC(:,1); |
| 148 |
if rem(ncdf,2) == 0, |
| 149 |
xAxU=[xG(:,1)' 2*xC(end,1)-xG(end,1)]'; |
| 150 |
else |
| 151 |
xAxU=S.XG(:,1); |
| 152 |
end |
| 153 |
else |
| 154 |
nx=varargin{4}; |
| 155 |
dxAx=max(max(xG(:)),max(xC(:)))-min(min(xG(:)),min(xC(:))); |
| 156 |
if dxAx > 360*(1-1/nx), dxAx=360; end |
| 157 |
dxAx=dxAx/nx; |
| 158 |
xAxU=min(xG(:))+[0:nx]'*dxAx; |
| 159 |
xAxC=(xAxU(1:nx)+xAxU(2:nx+1))/2; |
| 160 |
end |
| 161 |
|
| 162 |
%- volume: |
| 163 |
if rem(ncdf,2) == 1 | ktyp == 0, |
| 164 |
%fprintf(' rAc:'); fprintf(' %i',size(rAc)); fprintf('\n'); |
| 165 |
%fprintf('dims:'); fprintf(' %i',dims); fprintf('\n'); |
| 166 |
%fprintf(' dRf:'); fprintf(' %i',size(dRf)); fprintf('\n'); |
| 167 |
%vol=repmat(rAc,[1 1 dims(3)]).*hFacC; |
| 168 |
vol=reshape(rAc,[dims(1)*dims(2) 1])*dRf'; vol=reshape(vol,dims).*hFacC; |
| 169 |
end |
| 170 |
|
| 171 |
% create the structure |
| 172 |
|
| 173 |
if rem(ncdf,2) == 1 | ktyp == 0, |
| 174 |
G = struct('dims',dims, ... |
| 175 |
'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ... |
| 176 |
'xC',xC,'yC',yC,'xG',xG,'yG',yG,'rC',rC,'rF',rF, ... |
| 177 |
'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg,'dRc',dRc,'dRf',dRf, ... |
| 178 |
'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ... |
| 179 |
'csAngle',csAngle,'snAngle',snAngle, ... |
| 180 |
'hFacC',hFacC,'hFacW',hFacW,'hFacS',hFacS,'depth',depth,'vol',vol); |
| 181 |
elseif ktyp == 1, |
| 182 |
G = struct('dims',dims, ... |
| 183 |
'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ... |
| 184 |
'xC',xC,'yC',yC,'xG',xG,'yG',yG, ... |
| 185 |
'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg, ... |
| 186 |
'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ... |
| 187 |
'csAngle',csAngle,'snAngle',snAngle); |
| 188 |
else |
| 189 |
G = struct('dims',dims, ... |
| 190 |
'rC',rC,'rF',rF, ... |
| 191 |
'dRc',dRc,'dRf',dRf); |
| 192 |
end |
| 193 |
|
| 194 |
return |
Parent Directory
| 
Revision Log
| 
Revision Graph