Warning, /verification/tutorial_global_oce_latlon/diags_matlab/mit_loadgrid.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
051ee7f715 Jean*0001 function grid = mit_loadgrid(varargin);
0002 % function grid = mit_loadgrid(dname);
0003 % load the geometry of an arbitrary run with the mitgcm
0004
0005 % Aug 15, 2002: fixed a bug (?): yc' -> yc
0006 %
0007 if nargin == 0
0008 dname = '.';
0009 else
0010 dname = varargin{1};
0011 end
0012 % tracer time step
0013 deltattracer = mit_getparm(fullfile(dname,'data'),'deltaTtracer');
0014 if isempty(deltattracer)
0015 error('deltaTtracer is empty')
0016 end
0017 gravity = mit_getparm(fullfile(dname,'data'),'gravity');
0018 if isempty(gravity);
0019 disp('assuming gravity = 9.81')
0020 gravity = 9.81;
0021 end
0022 rhoNil = mit_getparm(fullfile(dname,'data'),'rhonil');
0023 if isempty(rhoNil);
0024 rhoNil = mit_getparm(fullfile(dname,'data'),'rhoNil');
0025 end
0026 if isempty(rhoNil);
0027 rhoNil = mit_getparm(fullfile(dname,'data'),'rhoConst');
0028 end
0029 if isempty(rhoNil);
0030 rhoNil = mit_getparm(fullfile(dname,'data'),'rhoconst');
0031 end
0032 if isempty(rhoNil);
0033 disp('assuming rhoNil = 1035')
0034 rhoNil = 1035;
0035 end
0036 % determine vertical coordinates
0037 br = mit_getparm(fullfile(dname,'data'),'buoyancyRelation');
0038 if isempty(br)
0039 disp('assuming buoyancyRelation = OCEANIC')
0040 br = 'OCEANIC';
0041 end
0042 if ~strcmpi(br,'OCEANIC');
0043 pfac = 1/rhoNil/gravity;
0044 else
0045 pfac = 1;
0046 end
0047 % create masks for plotting
0048 hfac=rdmds(fullfile(dname,'hFacC')); hfacc=change(hfac,'==',0,NaN);
0049 hfac=rdmds(fullfile(dname,'hFacW')); hfacw=change(hfac,'==',0,NaN);
0050 hfac=rdmds(fullfile(dname,'hFacS')); hfacs=change(hfac,'==',0,NaN);
0051 clear hfac;
0052 cmask=change(hfacc,'>',0,1);
0053 umask=change(hfacw,'>',0,1);
0054 vmask=change(hfacs,'>',0,1);
0055
0056 [nx ny nz] = size(cmask);
0057
0058 % find the index of the deepest wet tracer-cell
0059 klowc = sum(change(cmask,'==',NaN,0),3);
0060
0061 % create grid parameters
0062 dz = mit_getparm(fullfile(dname,'data'),'delR');
0063 if isempty(dz);
0064 disp('trying delZ')
0065 dz = mit_getparm(fullfile(dname,'data'),'delZ');
0066 if isempty(dz)
0067 error('vertical grid cannot be established')
0068 end
0069 end
0070 if size(dz,1) == 1;
0071 dz = dz';
0072 end
0073 if length(dz) ~= nz
0074 error('dz could not be retrieved correctly from the data file')
0075 end
0076 dz = dz*pfac;
0077
0078 if ~strcmp(br,'OCEANIC');
0079 dz = dz(end:-1:1);
0080 end
0081 zgpsi = [0;cumsum(dz)];
0082 zc = .5*(zgpsi(1:nz)+zgpsi(2:nz+1));
0083 zg=zgpsi(1:nz);
0084
0085 xc = rdmds(fullfile(dname,'XC'));
0086 yc = rdmds(fullfile(dname,'YC'));
0087 xg = rdmds(fullfile(dname,'XG'));
0088 yg = rdmds(fullfile(dname,'YG'));
0089 dxc = rdmds(fullfile(dname,'DXC'));
0090 dyc = rdmds(fullfile(dname,'DYC'));
0091 dxg = rdmds(fullfile(dname,'DXG'));
0092 dyg = rdmds(fullfile(dname,'DYG'));
0093
0094 lonc = xc(:,1);
0095 latc = yc(1,:)';
0096 long = xg(:,1);
0097 latg = yg(1,:)';
0098
0099 % depth
0100 depth = rdmds(fullfile(dname,'Depth'))*pfac;
0101
0102 % current directory
0103 if strcmp(dname,'.') | strcmp(dname,'./')
0104 [pathstr,dirname,ext]=fileparts(pwd);
0105 dirname = [dirname ext];
0106 else
0107 [pathstr,dirname,ext]=fileparts(fullfile(pwd,dname));
0108 dirname = [dirname ext];
0109 end
0110
0111 ius = findstr(dirname,'_');
0112 if ~isempty(ius)
0113 for k=1:length(ius)
0114 dirname = [dirname(1:ius(k)-1) '\' dirname(ius(k):end)];
0115 ius = ius+1;
0116 end
0117 end
0118 %
0119 if ~strcmp(br,'OCEANIC');
0120 hfacc = hfacc(:,:,end:-1:1);
0121 hfacs = hfacs(:,:,end:-1:1);
0122 hfacw = hfacw(:,:,end:-1:1);
0123 cmask = cmask(:,:,end:-1:1);
0124 umask = umask(:,:,end:-1:1);
0125 vmask = vmask(:,:,end:-1:1);
0126 klowc = nz - klowc + 1;
0127 klowc(find(klowc==nz+1)) = 0;
0128 end
0129
0130 % area and volume of grid cells has to be done after flipping the masks
0131 % vertically
0132 rac = rdmds(fullfile(dname,'RAC')).*cmask(:,:,1);
0133 ras = rdmds(fullfile(dname,'RAS')).*vmask(:,:,1);
0134 raw = rdmds(fullfile(dname,'RAW')).*umask(:,:,1);
0135 rac3d = repmat(rac,[1 1 nz]).*hfacc;
0136 volc = rac3d.*permute(repmat(dz,[1 nx ny]),[2 3 1]);
0137 % create the structure
0138 % Aug 15, 2002: fixed a bug (?): yc' -> yc
0139 grid = struct('dname',dirname, ...
0140 'deltattracer',deltattracer, ...
0141 'gravity',gravity','rhonil',rhoNil, ...
0142 'buoyancy',br,'pfac',pfac, ...
0143 'nx',nx,'ny',ny,'nz',nz, ...
0144 'xc',xc,'yc',yc,'xg',xg,'yg',yg, ...
0145 'lonc',lonc,'latc',latc,'long',long,'latg',latg, ...
0146 'dxc',dxc,'dyc',dyc,'dxg',dxg,'dyg',dyg, ...
0147 'rac',rac,'ras',ras,'raw',raw, ...
0148 'dz',dz,'zc',zc,'zg',zg,'zgpsi',zgpsi, ...
0149 'depth',depth,'hfacc',hfacc,'hfacs',hfacs,'hfacw',hfacw, ...
0150 'cmask',cmask,'umask',umask,'vmask',vmask,'volc',volc, ...
0151 'klowc',klowc);
0152
0153 return
