/[MITgcm]/MITgcm_contrib/gmaze_pv/eg_write_bin2cdf_latlongrid_subdomain.m
ViewVC logotype

Contents of /MITgcm_contrib/gmaze_pv/eg_write_bin2cdf_latlongrid_subdomain.m

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


Revision 1.3 - (show annotations) (download)
Wed Sep 19 15:37:38 2007 UTC (16 years, 7 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +0 -0 lines
General Update

1 % Script to extract and write in netCDF format a subdomain
2 % from the 1.8 global run
3 %
4 clear
5 global sla
6 pv_checkpath
7
8 % Load grid
9 GRID_125
10
11 % Load list of all outputs
12 otab = latlon8grid_outputs_table;
13
14 % Setup standard grid variables
15 lon_c = lon125;
16 lon_u = [lon125(1)-360+lon125(end) (lon125(2:end)+lon125(1:end-1))/2];
17 lat_c = lat125;
18 lat_v = [lat125(1)-(lat125(2)-lat125(1))/2 (lat125(1:end-1)+lat125(2:end))/2];
19 z_c = (cumsum(thk125)-thk125/2);
20 z_w = [0 cumsum(thk125(1:end-1))];
21
22
23 % Set subrange - Longitude given as degrees east
24 subdomain = 4;
25
26 switch subdomain
27 case 1
28 sub_name = 'western_north_atlantic';
29 lonmin = lon125(2209)-180;
30 lonmax = lon125(2497-1)-180;
31 latmin = lat125(1225);
32 latmax = lat125(1497-1);
33 depmin = min(z_w);
34 depmax = z_c(29);
35 m_proj('mercator','long',[270 365],'lat',[0 60]);
36 %clf;hold on;m_coast;m_grid;
37 LIMITS = [lonmin+180 lonmax+180 latmin latmax depmin depmax]
38 %m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','r','linewidth',2);
39 %title(sub_name);
40
41 case 3
42 sub_name = 'north_atlantic';
43 lonmin = lon125(2209)-180;
44 lonmax = lon125(2881-1)-180;
45 latmin = lat125(1157);
46 latmax = lat125(1565-1);
47 depmin = min(z_w);
48 depmax = z_c(29);
49 m_proj('mercator','long',[270 365],'lat',[0 60]);
50 clf;hold on;m_coast;m_grid;
51 LIMITS = [lonmin+180 lonmax+180 latmin latmax depmin depmax]
52 m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','k','linewidth',2);
53 title(sub_name);
54
55 case 4
56 sub_name = 'global';
57 lonmin = lon125(1)-180;
58 lonmax = lon125(2881-1)-180;
59 latmin = lat125(1);
60 latmax = lat125(2177-1);
61 depmin = min(z_w);
62 depmax = z_c(29); depmax = z_w(2);
63 m_proj('mercator','long',[0 360],'lat',[-90 90]);
64 clf;hold on;m_coast;m_grid;
65 LIMITS = [lonmin+180 lonmax+180 latmin latmax depmin depmax]
66 m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','k','linewidth',2);
67 title(sub_name);
68
69
70 case 10
71 sub_name = 'KE';
72 lonmin = lon125(961)-180;
73 lonmax = lon125(1601-1)-180;
74 latmin = lat125(1140);
75 latmax = lat125(1523-1);
76 depmin = min(z_w);
77 depmax = z_c(25);
78 m_proj('mercator','long',[0 360],'lat',[-90 90]);
79 %clf;hold on;m_coast;m_grid;
80 LIMITS = [lonmin+180 lonmax+180 latmin latmax depmin depmax]
81 %m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','k','linewidth',2);
82 %title(sub_name);
83
84
85 end
86
87 %refresh
88
89 % Path of the directory to find input binary files:
90 pathi = 'ecco2_cycle1_bin/';
91
92 % Path where the netcdf outputs will be stored:
93 patho = './ecco2_cycle1_netcdf/monthly/';
94 %patho = './ecco2_cycle1_netcdf/six_hourly/';
95
96 % Variables to analyse (index into otab):
97 wvar = [];
98 % 3D fields:
99 wvar = [wvar 34]; % THETA
100 wvar = [wvar 35]; % THETASQ
101 %wvar = [wvar 33]; % SALTanom
102 %wvar = [wvar 47]; % VVEL
103 %wvar = [wvar 31]; % RHOAnoma
104
105
106
107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 for i = 1 : length(wvar)
109 ifield = wvar(i);
110 fil = otab{ifield,1};
111 l = dir(strcat(pathi,fil,sla));
112 if ifield == 33,
113 l = dir(strcat(pathi,'SALT',sla));
114 end
115 if ifield == 35,
116 l = dir(strcat(pathi,'THETA',sla));
117 end
118 if ifield == 31,
119 l = dir(strcat(pathi,'RHO',sla));
120 end
121 it = 0;
122 clear ll
123 for il = 1 : size(l,1)
124 if ~l(il).isdir & findstr(l(il).name,strcat(fil,'.')) % is'it the file type we want ?
125 it = it + 1;
126 ll(it).name = l(il).name;
127 end %if
128 end %for il
129
130 if it ~= 0 % if we found any files to compute:
131
132 % Create the timetable:
133 for il = 1 : size(ll,2)
134 filinprog = ll(il).name;
135 stepnum=str2num(filinprog(findstr(filinprog,fil)+length(fil)+1:length(filinprog)- ...
136 length('.data')));
137 TIME(il,:) = dtecco2(stepnum,0);
138 end
139
140 % Translate files:
141 for il = 1 : size(ll,2)
142
143 filinprog = ll(il).name;
144 stepnum=str2num(filinprog(findstr(filinprog,fil)+length(fil)+1:length(filinprog)- ...
145 length('.data')));
146 ID = datenum(1992,1,1)+stepnum*300/60/60/24;
147 dte = datestr(ID,'yyyymmddHHMM');
148 disp(strcat(fil,'->',datestr(ID,'yyyy/mm/dd/HH:MM'),'== Recorded in ==>',TIME(il,:)));
149 dirout = strcat(patho,sla,TIME(il,:));
150
151 if 1 % Want to record ?
152 if ~exist(dirout,'dir')
153 mkdir(dirout);
154 end
155 pathname = strcat(pathi,fil,sla);
156 if ifield == 33,
157 pathname = strcat(pathi,'SALT',sla);
158 end
159 if ifield == 35,
160 pathname = strcat(pathi,'THETA',sla);
161 end
162 if ifield == 31,
163 pathname = strcat(pathi,'RHO',sla);
164 end
165 if ~exist(strcat(dirout,sla,sprintf('%s.%s.nc',otab{ifield,1},sub_name)),'file')
166 %if 1
167 latlon2ingrid_netcdf(pathname,strcat(dirout,sla),...
168 stepnum,otab{ifield,1},otab, ...
169 lon_c, lon_u, ...
170 lat_c, lat_v, ...
171 z_c, z_w, ...
172 sub_name, ...
173 lonmin,lonmax,latmin,latmax,depmin,depmax);
174 else
175 disp(strcat('##### Skip file (already done):',dirout,sla,...
176 sprintf('%s.%s.nc',otab{ifield,1},sub_name)))
177 end %if %file exist
178
179 end %if 1/0 want to record ?
180 % if il==1,break,end;
181
182 fclose('all');
183
184 end %for il
185
186 end %if it
187
188 end %for i

  ViewVC Help
Powered by ViewVC 1.1.22