/[MITgcm]/MITgcm/verification/global_ocean.90x40x15/diags_matlab/mit_loadgrid.m
ViewVC logotype

Annotation of /MITgcm/verification/global_ocean.90x40x15/diags_matlab/mit_loadgrid.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Sat Aug 12 20:27:17 2006 UTC (17 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
files moved to verification/tutorial_global_oce_latlon/diags_matlab

1 mlosch 1.1 function grid = mit_loadgrid(varargin);
2     % function grid = mit_loadgrid(dname);
3     % load the geometry of an arbitrary run with the mitgcm
4    
5     % Aug 15, 2002: fixed a bug (?): yc' -> yc
6    
7     if nargin == 0
8     dname = '.';
9     else
10     dname = varargin{1};
11     end
12     % tracer time step
13     deltattracer = mit_getparm(fullfile(dname,'data'),'deltaTtracer');
14     if isempty(deltattracer)
15     error('deltaTtracer is empty')
16     end
17     gravity = mit_getparm(fullfile(dname,'data'),'gravity');
18     if isempty(gravity);
19     disp('assuming gravity = 9.81')
20     gravity = 9.81;
21     end
22     rhoNil = mit_getparm(fullfile(dname,'data'),'rhonil');
23     if isempty(rhoNil);
24     rhoNil = mit_getparm(fullfile(dname,'data'),'rhoNil');
25     end
26     if isempty(rhoNil);
27     rhoNil = mit_getparm(fullfile(dname,'data'),'rhoConst');
28     end
29     if isempty(rhoNil);
30     rhoNil = mit_getparm(fullfile(dname,'data'),'rhoconst');
31     end
32     if isempty(rhoNil);
33     disp('assuming rhoNil = 1035')
34     rhoNil = 1035;
35     end
36     % determine vertical coordinates
37     br = mit_getparm(fullfile(dname,'data'),'buoyancyRelation');
38     if isempty(br)
39     disp('assuming buoyancyRelation = OCEANIC')
40     br = 'OCEANIC';
41     end
42     if ~strcmpi(br,'OCEANIC');
43     pfac = 1/rhoNil/gravity;
44     else
45     pfac = 1;
46     end
47     % create masks for plotting
48     hfac=rdmds(fullfile(dname,'hFacC')); hfacc=change(hfac,'==',0,NaN);
49     hfac=rdmds(fullfile(dname,'hFacW')); hfacw=change(hfac,'==',0,NaN);
50     hfac=rdmds(fullfile(dname,'hFacS')); hfacs=change(hfac,'==',0,NaN);
51     clear hfac;
52     cmask=change(hfacc,'>',0,1);
53     umask=change(hfacw,'>',0,1);
54     vmask=change(hfacs,'>',0,1);
55    
56     [nx ny nz] = size(cmask);
57    
58     % find the index of the deepest wet tracer-cell
59     klowc = sum(change(cmask,'==',NaN,0),3);
60    
61     % create grid parameters
62     dz = mit_getparm(fullfile(dname,'data'),'delR');
63     if isempty(dz);
64     disp('trying delZ')
65     dz = mit_getparm(fullfile(dname,'data'),'delZ');
66     if isempty(dz)
67     error('vertical grid cannot be established')
68     end
69     end
70     if size(dz,1) == 1;
71     dz = dz';
72     end
73     if length(dz) ~= nz
74     error('dz could not be retrieved correctly from the data file')
75     end
76     dz = dz*pfac;
77    
78     if ~strcmp(br,'OCEANIC');
79     dz = dz(end:-1:1);
80     end
81     zgpsi = [0;cumsum(dz)];
82     zc = .5*(zgpsi(1:nz)+zgpsi(2:nz+1));
83     zg=zgpsi(1:nz);
84    
85     xc = rdmds(fullfile(dname,'XC'));
86     yc = rdmds(fullfile(dname,'YC'));
87     xg = rdmds(fullfile(dname,'XG'));
88     yg = rdmds(fullfile(dname,'YG'));
89     dxc = rdmds(fullfile(dname,'DXC'));
90     dyc = rdmds(fullfile(dname,'DYC'));
91     dxg = rdmds(fullfile(dname,'DXG'));
92     dyg = rdmds(fullfile(dname,'DYG'));
93    
94     lonc = xc(:,1);
95     latc = yc(1,:)';
96     long = xg(:,1);
97     latg = yg(1,:)';
98    
99     % depth
100     depth = rdmds(fullfile(dname,'Depth'))*pfac;
101    
102     % current directory
103     if strcmp(dname,'.') | strcmp(dname,'./')
104     [pathstr,dirname,ext]=fileparts(pwd);
105     dirname = [dirname ext];
106     else
107     [pathstr,dirname,ext]=fileparts(fullfile(pwd,dname));
108     dirname = [dirname ext];
109     end
110    
111     ius = findstr(dirname,'_');
112     if ~isempty(ius)
113     for k=1:length(ius)
114     dirname = [dirname(1:ius(k)-1) '\' dirname(ius(k):end)];
115     ius = ius+1;
116     end
117     end
118     %
119     if ~strcmp(br,'OCEANIC');
120     hfacc = hfacc(:,:,end:-1:1);
121     hfacs = hfacs(:,:,end:-1:1);
122     hfacw = hfacw(:,:,end:-1:1);
123     cmask = cmask(:,:,end:-1:1);
124     umask = umask(:,:,end:-1:1);
125     vmask = vmask(:,:,end:-1:1);
126     klowc = nz - klowc + 1;
127     klowc(find(klowc==nz+1)) = 0;
128     end
129    
130     % area and volume of grid cells has to be done after flipping the masks
131     % vertically
132     rac = rdmds(fullfile(dname,'RAC')).*cmask(:,:,1);
133     ras = rdmds(fullfile(dname,'RAS')).*vmask(:,:,1);
134     raw = rdmds(fullfile(dname,'RAW')).*umask(:,:,1);
135     rac3d = repmat(rac,[1 1 nz]).*hfacc;
136     volc = rac3d.*permute(repmat(dz,[1 nx ny]),[2 3 1]);
137     % create the structure
138     % Aug 15, 2002: fixed a bug (?): yc' -> yc
139     grid = struct('dname',dirname, ...
140     'deltattracer',deltattracer, ...
141     'gravity',gravity','rhonil',rhoNil, ...
142     'buoyancy',br,'pfac',pfac, ...
143     'nx',nx,'ny',ny,'nz',nz, ...
144     'xc',xc,'yc',yc,'xg',xg,'yg',yg, ...
145     'lonc',lonc,'latc',latc,'long',long,'latg',latg, ...
146     'dxc',dxc,'dyc',dyc,'dxg',dxg,'dyg',dyg, ...
147     'rac',rac,'ras',ras,'raw',raw, ...
148     'dz',dz,'zc',zc,'zg',zg,'zgpsi',zgpsi, ...
149     'depth',depth,'hfacc',hfacc,'hfacs',hfacs,'hfacw',hfacw, ...
150     'cmask',cmask,'umask',umask,'vmask',vmask,'volc',volc, ...
151     'klowc',klowc);
152    
153     return

  ViewVC Help
Powered by ViewVC 1.1.22