Warning, /utils/matlab/load_grid.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit e6feef68 on 2025-11-18 19:40:16 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
e6feef680a Ivan*0030 % kinp: = 0 : binary file ; = 1 : old MNC file ; = 2 : new MNC file
ecaf506031 Jean*0031
0032 if ncdf > 1, TimeT0=clock; end
0033
0034 if rem(ncdf,2) == 0,
0035 %- load MDSIO grid files :
e6feef680a Ivan*0036 kinp=0;
ecaf506031 Jean*0037 if ktyp < 2,
0038 xC=rdmds(fullfile(rDir,'XC'));
0039 yC=rdmds(fullfile(rDir,'YC'));
0040 xG=rdmds(fullfile(rDir,'XG'));
0041 yG=rdmds(fullfile(rDir,'YG'));
0042
0043 dXc=rdmds(fullfile(rDir,'DXC'));
0044 dYc=rdmds(fullfile(rDir,'DYC'));
0045 dXg=rdmds(fullfile(rDir,'DXG'));
0046 dYg=rdmds(fullfile(rDir,'DYG'));
0047
0048 rAc=rdmds(fullfile(rDir,'RAC'));
0049 rAw=rdmds(fullfile(rDir,'RAW'));
0050 rAs=rdmds(fullfile(rDir,'RAS'));
0051 rAz=rdmds(fullfile(rDir,'RAZ'));
0052
0053 %- load grid orientation relative to West-East / South-North dir:
0054 if isempty(dir(fullfile(rDir,'AngleCS.*'))),
0055 fprintf(' no grid orientation (Cos & Sin) file to load ');
0056 csAngle= ones(size(rAc));
0057 snAngle=zeros(size(rAc));
0058 fprintf(' => set COS=1,SIN=0\n');
0059 else
0060 csAngle=rdmds(fullfile(rDir,'AngleCS'));
0061 snAngle=rdmds(fullfile(rDir,'AngleSN'));
0062 end
0063 dims = size(rAc);
0064 end
0065
0066 if rem(ktyp,2) == 0,
0067 rC=rdmds(fullfile(rDir,'RC')); rC=squeeze(rC);
0068 rF=rdmds(fullfile(rDir,'RF')); rF=squeeze(rF);
0069 dRc=rdmds(fullfile(rDir,'DRC')); dRc=squeeze(dRc);
0070 dRf=rdmds(fullfile(rDir,'DRF')); dRf=squeeze(dRf);
6b4f212f67 Jean*0071 depth=rdmds(fullfile(rDir,'Depth'));
0072 dims = [size(depth) length(rC)];
ecaf506031 Jean*0073 end
0074
0075 if ktyp == 0,
0076 %- define domain:
0077 hFacC=rdmds(fullfile(rDir,'hFacC'));
0078 hFacW=rdmds(fullfile(rDir,'hFacW'));
0079 hFacS=rdmds(fullfile(rDir,'hFacS'));
0080 dims = size(hFacC);
0081 end
0082
0083 if ncdf > 1, TimeT1=clock; end
0084
0085 else
0086 %- load NetCDF grid files :
0087 S=rdmnc(fullfile(rDir,'grid.*.nc'));
0088 if ncdf > 1, TimeT1=clock; end
0089
6b4f212f67 Jean*0090 %--
e6feef680a Ivan*0091 if isfield(S,'XC'),
0092 kinp=1; % old MNC grid names until c69h
0093 xC=S.XC;
0094 yC=S.YC;
0095 xG=S.XG(1:end-1,1:end-1);
0096 yG=S.YG(1:end-1,1:end-1);
0097 rC=S.RC';
0098 rF=S.RF';
0099 rAc=S.rA;
0100 if isfield(S,'AngleCS') & isfield(S,'AngleSN'),
0101 csAngle=S.AngleCS;
0102 snAngle=S.AngleSN;
0103 else
0104 fprintf(' no grid orientation (Cos & Sin) in grid file ');
0105 csAngle= ones(size(rAc));
0106 snAngle=zeros(size(rAc));
0107 fprintf(' => set COS=1,SIN=0\n');
0108 end
0109 hFacC=S.HFacC;
0110 hFacW=S.HFacW(1:end-1,:,:);
0111 hFacS=S.HFacS(:,1:end-1,:);
0112 else
0113 kinp=2; % current MNC grid names >=c69h
0114 xC=S.xC;
0115 yC=S.yC;
0116 xG=S.xG(1:end-1,1:end-1);
0117 yG=S.yG(1:end-1,1:end-1);
0118 rC=S.rC';
0119 rF=S.rF';
0120 rAc=S.rAc;
0121 if isfield(S,'angleCS') & isfield(S,'angleSN'),
0122 csAngle=S.angleCS;
0123 snAngle=S.angleSN;
0124 else
0125 fprintf(' no grid orientation (Cos & Sin) in grid file ');
0126 csAngle= ones(size(rAc));
0127 snAngle=zeros(size(rAc));
0128 fprintf(' => set COS=1,SIN=0\n');
0129 end
0130 hFacC=S.hFacC;
0131 hFacW=S.hFacW(1:end-1,:,:);
0132 hFacS=S.hFacS(:,1:end-1,:);
0133 end
ecaf506031 Jean*0134 dXc=S.dxC(1:end-1,:);
0135 dYc=S.dyC(:,1:end-1);
0136 dXg=S.dxG(:,1:end-1);
0137 dYg=S.dyG(1:end-1,:);
3d4f3e5a15 Jean*0138 dRc=S.drC';
0139 dRf=S.drF';
ecaf506031 Jean*0140 rAw=S.rAw(1:end-1,:);
0141 rAs=S.rAs(:,1:end-1);
0142 rAz=S.rAz(1:end-1,1:end-1);
0143 depth=S.Depth;
0144 dims = size(hFacC);
e6feef680a Ivan*0145
ecaf506031 Jean*0146 end
0147
e6feef680a Ivan*0148 if ncdf > 1, fprintf(' (took %6.4f s)',etime(TimeT1,TimeT0)); end
ecaf506031 Jean*0149
0150 if length(dims) == 1, dims(2)=1; end
0151 if length(dims) == 2, dims(3)=1; end
0152
6b4f212f67 Jean*0153 if rem(ncdf,2) == 1 | ktyp ~= 2,
ecaf506031 Jean*0154 %- 1-D axis:
6b4f212f67 Jean*0155 if nargin < 3,
ecaf506031 Jean*0156 %- yAxis
6b4f212f67 Jean*0157 ny=dims(2);
0158 yAxC=yC(1,:)';
e6feef680a Ivan*0159 if kinp == 0,
6b4f212f67 Jean*0160 yAxV=[yG(1,:) 2*yC(1,end)-yG(1,end)]';
e6feef680a Ivan*0161 elseif kinp == 1,
6b4f212f67 Jean*0162 yAxV=S.YG(1,:)';
e6feef680a Ivan*0163 else
0164 yAxV=S.yG(1,:)';
6b4f212f67 Jean*0165 end
ecaf506031 Jean*0166 else
6b4f212f67 Jean*0167 ny=varargin{3};
0168 dyAx=(max(yG(:))-min(yG(:)))/ny;
0169 yAxV=min(yG(:))+[0:ny]'*dyAx;
0170 yAxC=(yAxV(1:ny)+yAxV(2:ny+1))/2;
ecaf506031 Jean*0171 end
6b4f212f67 Jean*0172 if nargin < 4,
0173 %- xAxis
0174 nx=dims(1);
0175 xAxC=xC(:,1);
e6feef680a Ivan*0176 if kinp == 0,
6b4f212f67 Jean*0177 xAxU=[xG(:,1)' 2*xC(end,1)-xG(end,1)]';
e6feef680a Ivan*0178 elseif kinp == 1,
0179 xAxU=S.XG(:,1);
6b4f212f67 Jean*0180 else
e6feef680a Ivan*0181 xAxU=S.xG(:,1);
6b4f212f67 Jean*0182 end
ecaf506031 Jean*0183 else
6b4f212f67 Jean*0184 nx=varargin{4};
0185 dxAx=max(max(xG(:)),max(xC(:)))-min(min(xG(:)),min(xC(:)));
0186 if dxAx > 360*(1-1/nx), dxAx=360; end
0187 dxAx=dxAx/nx;
0188 xAxU=min(xG(:))+[0:nx]'*dxAx;
0189 xAxC=(xAxU(1:nx)+xAxU(2:nx+1))/2;
ecaf506031 Jean*0190 end
0191
0192 %- volume:
6b4f212f67 Jean*0193 if rem(ncdf,2) == 1 | ktyp == 0,
0194 vol=reshape(rAc,[dims(1)*dims(2) 1])*dRf'; vol=reshape(vol,dims).*hFacC;
0195 end
ecaf506031 Jean*0196 end
0197
6b4f212f67 Jean*0198 %- clear space:
0199 if rem(ncdf,2) == 1, clear S ; end
0200
ecaf506031 Jean*0201 % create the structure
0202
0203 if rem(ncdf,2) == 1 | ktyp == 0,
0204 G = struct('dims',dims, ...
0205 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
0206 'xC',xC,'yC',yC,'xG',xG,'yG',yG,'rC',rC,'rF',rF, ...
0207 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg,'dRc',dRc,'dRf',dRf, ...
0208 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
0209 'csAngle',csAngle,'snAngle',snAngle, ...
0210 'hFacC',hFacC,'hFacW',hFacW,'hFacS',hFacS,'depth',depth,'vol',vol);
0211 elseif ktyp == 1,
0212 G = struct('dims',dims, ...
0213 'nx',nx,'ny',ny,'xAxC',xAxC,'yAxC',yAxC,'xAxU',xAxU,'yAxV',yAxV, ...
0214 'xC',xC,'yC',yC,'xG',xG,'yG',yG, ...
0215 'dXc',dXc,'dYc',dYc,'dXg',dXg,'dYg',dYg, ...
0216 'rAc',rAc,'rAw',rAw,'rAs',rAs,'rAz',rAz, ...
0217 'csAngle',csAngle,'snAngle',snAngle);
0218 else
0219 G = struct('dims',dims, ...
0220 'rC',rC,'rF',rF, ...
6b4f212f67 Jean*0221 'dRc',dRc,'dRf',dRf,'depth',depth);
ecaf506031 Jean*0222 end
0223
6b4f212f67 Jean*0224 fprintf(' and leaving\n');
0225 %return
0226 end
