Warning, /utils/matlab/load_grid.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 6b4f212f on 2024-06-05 16:55:17 UTC
6b4f212f67 Jean*0001 function G = load_grid(varargin)
ecaf506031 Jean*0002 % function G = load_grid(dirName,option,nyAxis,nxAxis);
0003 %
0004 % load the MITgcm output grid-files into one structure array
0005 % arguments:
0006 % dirName = directory where grid-files are
0007 % optional arguments:
0008 % ncdf = last digit of integer "option":
0009 % ncdf = 0,2 : read binary (MDSIO) files using rdmds (default: ncdf=0)
0010 % ncdf = 1,3 : read NetCDF (MNC) files using rdmnc
0011 % ncdf = 2,3 : as 0,1 + print elapsed-time for loading files
0012 % option >= 10 : only read in Horiz.Grid spacing (without hFac)
0013 % option >= 20 : only read in Verti.Grid spacing (without hFac)
6b4f212f67 Jean*0014 % nyAxis = number of spacing in 1.D y-axis yAxC (default: nyAxis = grid-Ny)
ecaf506031 Jean*0015 % nxAxis = number of spacing in 1.D x-axis xAxC (default: nxAxis = grid-Nx)
0016
6b4f212f67 Jean*0017 fprintf('entering fct "load_grid" ...');
ecaf506031 Jean*0018 if nargin == 0
0019 rDir = '.';
0020 else
0021 rDir = varargin{1};
0022 end
0023 if nargin < 2,
0024 ktyp=0;
0025 ncdf=0;
0026 else
0027 ncdf= varargin{2};
0028 ktyp=floor(ncdf/10); ncdf=rem(ncdf,10);
0029 end
0030
0031 if ncdf > 1, TimeT0=clock; end
0032
0033 if rem(ncdf,2) == 0,
0034 %- load MDSIO grid files :
0035 if ktyp < 2,
0036 xC=rdmds(fullfile(rDir,'XC'));
0037 yC=rdmds(fullfile(rDir,'YC'));
0038 xG=rdmds(fullfile(rDir,'XG'));
0039 yG=rdmds(fullfile(rDir,'YG'));
0040
0041 dXc=rdmds(fullfile(rDir,'DXC'));
0042 dYc=rdmds(fullfile(rDir,'DYC'));
0043 dXg=rdmds(fullfile(rDir,'DXG'));
0044 dYg=rdmds(fullfile(rDir,'DYG'));
0045
0046 rAc=rdmds(fullfile(rDir,'RAC'));
0047 rAw=rdmds(fullfile(rDir,'RAW'));
0048 rAs=rdmds(fullfile(rDir,'RAS'));
0049 rAz=rdmds(fullfile(rDir,'RAZ'));
0050
0051 %- load grid orientation relative to West-East / South-North dir:
0052 if isempty(dir(fullfile(rDir,'AngleCS.*'))),
0053 fprintf(' no grid orientation (Cos & Sin) file to load ');
0054 csAngle= ones(size(rAc));
0055 snAngle=zeros(size(rAc));
0056 fprintf(' => set COS=1,SIN=0\n');
0057 else
0058 csAngle=rdmds(fullfile(rDir,'AngleCS'));
0059 snAngle=rdmds(fullfile(rDir,'AngleSN'));
0060 end
0061 dims = size(rAc);
0062 end
0063
0064 if rem(ktyp,2) == 0,
0065 rC=rdmds(fullfile(rDir,'RC')); rC=squeeze(rC);
0066 rF=rdmds(fullfile(rDir,'RF')); rF=squeeze(rF);
0067 dRc=rdmds(fullfile(rDir,'DRC')); dRc=squeeze(dRc);
0068 dRf=rdmds(fullfile(rDir,'DRF')); dRf=squeeze(dRf);
6b4f212f67 Jean*0069 depth=rdmds(fullfile(rDir,'Depth'));
0070 dims = [size(depth) length(rC)];
ecaf506031 Jean*0071 end
0072
0073 if ktyp == 0,
0074 %- define domain:
0075 hFacC=rdmds(fullfile(rDir,'hFacC'));
0076 hFacW=rdmds(fullfile(rDir,'hFacW'));
0077 hFacS=rdmds(fullfile(rDir,'hFacS'));
0078 dims = size(hFacC);
0079 end
0080
0081 if ncdf > 1, TimeT1=clock; end
0082
0083 else
0084 %- load NetCDF grid files :
0085 %S=rdmnc([rDir,'grid.*.nc']);
0086 S=rdmnc(fullfile(rDir,'grid.*.nc'));
0087 if ncdf > 1, TimeT1=clock; end
0088
6b4f212f67 Jean*0089 %--
ecaf506031 Jean*0090 xC=S.XC;
0091 yC=S.YC;
0092 xG=S.XG(1:end-1,1:end-1);
0093 yG=S.YG(1:end-1,1:end-1);
3d4f3e5a15 Jean*0094 rC=S.RC';
0095 rF=S.RF';
ecaf506031 Jean*0096 dXc=S.dxC(1:end-1,:);
0097 dYc=S.dyC(:,1:end-1);
0098 dXg=S.dxG(:,1:end-1);
0099 dYg=S.dyG(1:end-1,:);
3d4f3e5a15 Jean*0100 dRc=S.drC';
0101 dRf=S.drF';
ecaf506031 Jean*0102 rAc=S.rA;
0103 rAw=S.rAw(1:end-1,:);
0104 rAs=S.rAs(:,1:end-1);
0105 rAz=S.rAz(1:end-1,1:end-1);
6b4f212f67 Jean*0106 if isfield(S,'AngleCS') & isfield(S,'AngleSN'),
ecaf506031 Jean*0107 csAngle=S.AngleCS;
0108 snAngle=S.AngleSN;
0109 else
0110 fprintf(' no grid orientation (Cos & Sin) in grid file ');
0111 csAngle= ones(size(rAc));
0112 snAngle=zeros(size(rAc));
0113 fprintf(' => set COS=1,SIN=0\n');
0114 end
0115 hFacC=S.HFacC;
0116 hFacW=S.HFacW(1:end-1,:,:);
0117 hFacS=S.HFacS(:,1:end-1,:);
0118 depth=S.Depth;
0119 dims = size(hFacC);
0120 end
0121
0122 if ncdf > 1, fprintf(' (took %6.4f s)\n',etime(TimeT1,TimeT0)); end
0123
0124 if length(dims) == 1, dims(2)=1; end
0125 if length(dims) == 2, dims(3)=1; end
0126
6b4f212f67 Jean*0127 if rem(ncdf,2) == 1 | ktyp ~= 2,
ecaf506031 Jean*0128 %- 1-D axis:
6b4f212f67 Jean*0129 if nargin < 3,
ecaf506031 Jean*0130 %- yAxis
6b4f212f67 Jean*0131 ny=dims(2);
0132 yAxC=yC(1,:)';
0133 if rem(ncdf,2) == 0,
0134 yAxV=[yG(1,:) 2*yC(1,end)-yG(1,end)]';
0135 else
0136 yAxV=S.YG(1,:)';
0137 end
ecaf506031 Jean*0138 else
6b4f212f67 Jean*0139 ny=varargin{3};
0140 dyAx=(max(yG(:))-min(yG(:)))/ny;
0141 yAxV=min(yG(:))+[0:ny]'*dyAx;
0142 yAxC=(yAxV(1:ny)+yAxV(2:ny+1))/2;
ecaf506031 Jean*0143 end
6b4f212f67 Jean*0144 if nargin < 4,
0145 %- xAxis
0146 nx=dims(1);
0147 xAxC=xC(:,1);
0148 if rem(ncdf,2) == 0,
0149 xAxU=[xG(:,1)' 2*xC(end,1)-xG(end,1)]';
0150 else
0151 xAxU=S.XG(:,1);
0152 end
ecaf506031 Jean*0153 else
6b4f212f67 Jean*0154 nx=varargin{4};
0155 dxAx=max(max(xG(:)),max(xC(:)))-min(min(xG(:)),min(xC(:)));
0156 if dxAx > 360*(1-1/nx), dxAx=360; end
0157 dxAx=dxAx/nx;
0158 xAxU=min(xG(:))+[0:nx]'*dxAx;
0159 xAxC=(xAxU(1:nx)+xAxU(2:nx+1))/2;
ecaf506031 Jean*0160 end
0161
0162 %- volume:
6b4f212f67 Jean*0163 if rem(ncdf,2) == 1 | ktyp == 0,
0164 vol=reshape(rAc,[dims(1)*dims(2) 1])*dRf'; vol=reshape(vol,dims).*hFacC;
0165 end
ecaf506031 Jean*0166 end
0167
6b4f212f67 Jean*0168 %- clear space:
0169 if rem(ncdf,2) == 1, clear S ; end
0170
ecaf506031 Jean*0171 % create the structure
0172
0173 if rem(ncdf,2) == 1 | ktyp == 0,
0174 G = struct('dims',dims, ...
0175 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
0176 'xC',xC,'yC',yC,'xG',xG,'yG',yG,'rC',rC,'rF',rF, ...
0177 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg,'dRc',dRc,'dRf',dRf, ...
0178 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
0179 'csAngle',csAngle,'snAngle',snAngle, ...
0180 'hFacC',hFacC,'hFacW',hFacW,'hFacS',hFacS,'depth',depth,'vol',vol);
0181 elseif ktyp == 1,
0182 G = struct('dims',dims, ...
0183 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
0184 'xC',xC,'yC',yC,'xG',xG,'yG',yG, ...
0185 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg, ...
0186 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
0187 'csAngle',csAngle,'snAngle',snAngle);
0188 else
0189 G = struct('dims',dims, ...
0190 'rC',rC,'rF',rF, ...
6b4f212f67 Jean*0191 'dRc',dRc,'dRf',dRf,'depth',depth);
ecaf506031 Jean*0192 end
0193
6b4f212f67 Jean*0194 fprintf(' and leaving\n');
0195 %return
0196 end
