MITgcm_contrib/gmaze_pv/eg_write_bin2cdf_csgrid_subdomain.m

% Script to extract a subdomain from a CS510 simulation
% and write in netCDF format on a regular lat/lon grid (1/4)
%
clear
global sla
pv_checkpath


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global setup:
% Restrict global domain to:
subdomain = 3; % North Atlantic

% Path to find input binary Cube sphere files:
pathi = './bin_cube49';

% Path where the netcdf outputs will be stored:
patho = './ecco2_cube49_netcdf';
patho = strcat(patho,sla,'monthly');
%patho = strcat(patho,sla,'six_hourly');

% Time step (for date conversion):
dt = 1200;

% Variables to analyse (index into otab):
otab = cs510grid_outputs_table; % from the 1/8 latlon definition
wvar = [];
dimen = 3;
switch dimen
  case 3 % 3D fields:
    %wvar = [wvar 34]; % THETA
    %wvar = [wvar 31]; % RHOAnoma
    %wvar = [wvar 33]; % SALTanom
  case 2 % 2D fields:
    wvar = [wvar 23]; % TFLUX
    %wvar = [wvar 20]; % SST
    %wvar = [wvar 19]; % SSS
end %switch number of dimensions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pre-process
% Get the grid:
path_grid = './grid';
XC = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'XC')),1,'float32');
YC = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'YC')),1,'float32');
XG = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'XG')),1,'float32');
YG = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'YG')),1,'float32');
GRID_125;
ZC = - [0 cumsum(thk125(1:length(thk125)-1))];
clear dpt125 lat125 lon125 thk125

% How to move to a lat/lon grid:
% CS510 is about 22km average resolution, ie: 1/4 degree
XI =   -180 : 1/4 : 180; 
YI = -90 : 1/4 : 90;
ZI = ZC;
if ~exist('CS510_to_LATLON025.mat','file')
   del = cube2latlon_preprocess(XC,YC,XI,YI);
   save('CS510_to_LATLON025.mat','XI','YI','XC','YC','del','-v6');
else
   load('CS510_to_LATLON025.mat')
end

% Set subrange - Longitude given as degrees east 
% (exact values come from the 1/8 lat-lon grid)
switch subdomain
  case 3
   sub_name = 'north_atlantic';
   lonmin = 276.0625;
   lonmax = 359.9375;
   latmin = 12.0975;
   latmax = 53.2011;
   depmin = 1;    % !!! indices
   depmax = 29;   % !!! indices
   if dimen == 3, depmax = 29,
   else, depmax = 1;end
   LIMITS = [lonmin lonmax latmin latmax depmin depmax]
   if 0
     m_proj('mercator','long',[270 365],'lat',[0 60]);
     clf;hold on;m_coast;m_grid;
     m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','k','linewidth',2);
     title(sub_name);
   end %if 1/0
end

% Get subdomain horizontal axis:
xi = XI(max(find(XI<=LIMITS(1))):max(find(XI<=LIMITS(2))));
yi = YI(max(find(YI<=LIMITS(3))):max(find(YI<=LIMITS(4))));
zi = ZI(LIMITS(5):LIMITS(6));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over variables to read
for i = 1 : length(wvar)
 ifield = wvar(i);
 fil = otab{ifield,1};
 
%                                                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Get info over the time loop
 % Get the file list:
 fild = fil; % Insert condition here for special files:
 if ifield == 23 & findstr(patho,'six')
    fild = 'surForcT';
    fil  = 'surForcT';
 end   
 l = dir(strcat(pathi,sla,fild));
 it = 0;
 clear ll
 for il = 1 : size(l,1)
   if ~l(il).isdir & findstr(l(il).name,strcat(fil,'.')) % is'it the file type we want ?
     it = it + 1;
     ll(it).name = l(il).name;
   end %if
 end %for il
 % Create the timetable:
 for il = 1 : size(ll,2)
   filin = ll(il).name;
   % Now extract the stepnum from : %s.%10.10d.data
   ic = findstr(filin,fil)+length(fil)+1; i = 1; clear stepnum
   while filin(ic) ~= '.'
      stepnum(i) = filin(ic); i = i + 1; ic = ic + 1;
   end
   ID = str2num(stepnum);
   TIME(il,:) = datestr(datenum(1992,1,1)+ID*dt/60/60/24,'yyyymmddHHMM');
 end
 nt = size(TIME,1);
 
%                                                 %%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over time
 for it = 1 : nt
   snapshot = TIME(it,:);
   ID = 60*60*24/dt*( datenum(snapshot,'yyyymmddHHMM') - datenum(1992,1,1) );
   filin = ll(it).name;
   disp('')
   disp(strcat('Processing: ',fil,'//',snapshot))
   dirout = strcat(patho,sla,TIME(it,:),sla);
   filout = sprintf('%s.%s.nc',otab{ifield,1},sub_name);
   
   if ~exist(strcat(dirout,sla,filout),'file') % File already exists ?
   
%%%% READ THE FILE
   switch otab{ifield,6}
    case 4, flt = 'float32';
    case 8, flt = 'float64';
   end
   t0 = clock;
   if findstr(filin,'.gz') % Gunzip file, special care !     
      disp('|----> Find a file with gz extension, work on uncompressed file ...')
      
      % 1: copy the filename with it path into gunzip_1_file.txt
      fid1 = fopen('gunzip_1_file.txt','w');
      fprintf(fid1,'%s',strcat(pathi,sla,fild,sla,filin));fclose(fid1);

      % 2: uncompress the file into a temporary folder:
      disp('|------> uncompressing the file ...')
      ! ./gunzip_1_file.bat
      disp(strcat('|--------: ',num2str(etime(t0,clock))))
      
      % 3: Read the uncompress file:
      disp('|--> reading it ...')
      C = readrec_cs510(strcat('gunzip_1_file',sla,'tempo.data'),LIMITS(6),flt);
      disp(strcat('|----: ',num2str(etime(t0,clock))))
      
      % 4: Suppress it
      ! \rm ./gunzip_1_file/tempo.data
      
   else % Simply read the file:
      disp('|--> reading it ...')
      C = readrec_cs510(strcat(pathi,sla,fild,sla,filin),LIMITS(6),flt);
      disp(strcat('|----: ',num2str(etime(t0,clock))))
   end
   
%%%% RESTRICT TO SUBDOMAIN
   disp('|--> get subdomain ...')
   % Restrict vertical to subdomain:
   if LIMITS(5) ~= 1
      disp('|----> vertical ...');
      C = C(:,:,LIMITS(5):end);
   end
   % Clean the field:
   C(find(C==0)) = NaN; 
   % Move the field into lat/lon grid:
   disp('|----> Move to lat/lon grid ...');
   C = cube2latlon_fast(del,C);
   % And then restrict horizontal to subdomain: 
   disp('|----> horizontal ...');  
   C = C(max(find(XI<=LIMITS(1))):max(find(XI<=LIMITS(2))),...
         max(find(YI<=LIMITS(3))):max(find(YI<=LIMITS(4))),:);
   
   
%%%% RECORD
   disp('|--> record netcdf file ...')
   fid1 = fopen('inprogress.txt','w');
   fprintf(fid1,'%s',strcat(dirout,sla,filout));fclose(fid1);
   
   if 1 % Realy want to record ?
     
   if ~exist(dirout,'dir')
       mkdir(dirout);
   end
   
   
     nc = netcdf('inprogress.nc','clobber');

     nc('X') = length(xi);
     nc('Y') = length(yi);
     nc('Z') = length(zi);

     nc{'X'}='X';
     nc{'Y'}='Y';
     nc{'Z'}='Z';

     nc{'X'}.uniquename='X';
     nc{'X'}.long_name='longitude';
     nc{'X'}.gridtype=ncint(0);
     nc{'X'}.units='degrees_east';
     nc{'X'}(:) = xi;

     nc{'Y'}.uniquename='Y';
     nc{'Y'}.long_name='latitude';
     nc{'Y'}.gridtype=ncint(0);
     nc{'Y'}.units='degrees_north';
     nc{'Y'}(:) = yi;

     nc{'Z'}.uniquename='Z';
     nc{'Z'}.long_name='depth';
     nc{'Z'}.gridtype=ncint(0);
     nc{'Z'}.units='m';
     nc{'Z'}(:) = zi;

     ncid = fil;
     nc{ncid}={'Z' 'Y' 'X'};
     nc{ncid}.missing_value = ncdouble(NaN);
     nc{ncid}.FillValue_ = ncdouble(0.0);
     nc{ncid}(:,:,:) = permute(C,[3 2 1]);
     nc{ncid}.units = otab{ifield,5};

     close(nc);
     ! ./inprogress.bat
     
   end %if 1/0 want to record ?
   disp(strcat('|--: ',num2str(etime(t0,clock))))
   
   else
     disp(strcat('|--> Skip file (already done):',dirout,sla,filout))
   end %if %file exist
 
 end %for it
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END Loop over time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
end %if it
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END Loop over variables to read
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1	gmaze	1.1	% Script to extract a subdomain from a CS510 simulation
2			% and write in netCDF format on a regular lat/lon grid (1/4)
3			%
4			clear
5			global sla
6			pv_checkpath
7
8
9			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Global setup:
10			% Restrict global domain to:
11			subdomain = 3; % North Atlantic
12
13			% Path to find input binary Cube sphere files:
14			pathi = './bin_cube49';
15
16			% Path where the netcdf outputs will be stored:
17			patho = './ecco2_cube49_netcdf';
18			patho = strcat(patho,sla,'monthly');
19			%patho = strcat(patho,sla,'six_hourly');
20
21			% Time step (for date conversion):
22			dt = 1200;
23
24			% Variables to analyse (index into otab):
25			otab = cs510grid_outputs_table; % from the 1/8 latlon definition
26			wvar = [];
27			dimen = 3;
28			switch dimen
29			case 3 % 3D fields:
30			%wvar = [wvar 34]; % THETA
31			%wvar = [wvar 31]; % RHOAnoma
32			%wvar = [wvar 33]; % SALTanom
33			case 2 % 2D fields:
34			wvar = [wvar 23]; % TFLUX
35			%wvar = [wvar 20]; % SST
36			%wvar = [wvar 19]; % SSS
37			end %switch number of dimensions
38
39			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Pre-process
40			% Get the grid:
41			path_grid = './grid';
42			XC = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'XC')),1,'float32');
43			YC = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'YC')),1,'float32');
44			XG = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'XG')),1,'float32');
45			YG = readrec_cs510(sprintf('%s.data',strcat(path_grid,sla,'YG')),1,'float32');
46			GRID_125;
47			ZC = - [0 cumsum(thk125(1:length(thk125)-1))];
48			clear dpt125 lat125 lon125 thk125
49
50			% How to move to a lat/lon grid:
51			% CS510 is about 22km average resolution, ie: 1/4 degree
52			XI = -180 : 1/4 : 180;
53			YI = -90 : 1/4 : 90;
54			ZI = ZC;
55			if ~exist('CS510_to_LATLON025.mat','file')
56			del = cube2latlon_preprocess(XC,YC,XI,YI);
57			save('CS510_to_LATLON025.mat','XI','YI','XC','YC','del','-v6');
58			else
59			load('CS510_to_LATLON025.mat')
60			end
61
62			% Set subrange - Longitude given as degrees east
63			% (exact values come from the 1/8 lat-lon grid)
64			switch subdomain
65			case 3
66			sub_name = 'north_atlantic';
67			lonmin = 276.0625;
68			lonmax = 359.9375;
69			latmin = 12.0975;
70			latmax = 53.2011;
71			depmin = 1; % !!! indices
72			depmax = 29; % !!! indices
73			if dimen == 3, depmax = 29,
74			else, depmax = 1;end
75			LIMITS = [lonmin lonmax latmin latmax depmin depmax]
76			if 0
77			m_proj('mercator','long',[270 365],'lat',[0 60]);
78			clf;hold on;m_coast;m_grid;
79			m_line(LIMITS([1 2 2 1 1]),LIMITS([3 3 4 4 3]),'color','k','linewidth',2);
80			title(sub_name);
81			end %if 1/0
82			end
83
84			% Get subdomain horizontal axis:
85			xi = XI(max(find(XI<=LIMITS(1))):max(find(XI<=LIMITS(2))));
86			yi = YI(max(find(YI<=LIMITS(3))):max(find(YI<=LIMITS(4))));
87			zi = ZI(LIMITS(5):LIMITS(6));
88
89			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over variables to read
91			for i = 1 : length(wvar)
92			ifield = wvar(i);
93			fil = otab{ifield,1};
94
95			% %%%%%%%%%%%%%%%%%%%%%%%%%%%%
96			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Get info over the time loop
97			% Get the file list:
98			fild = fil; % Insert condition here for special files:
99			if ifield == 23 & findstr(patho,'six')
100			fild = 'surForcT';
101			fil = 'surForcT';
102			end
103			l = dir(strcat(pathi,sla,fild));
104			it = 0;
105			clear ll
106			for il = 1 : size(l,1)
107			if ~l(il).isdir & findstr(l(il).name,strcat(fil,'.')) % is'it the file type we want ?
108			it = it + 1;
109			ll(it).name = l(il).name;
110			end %if
111			end %for il
112			% Create the timetable:
113			for il = 1 : size(ll,2)
114			filin = ll(il).name;
115			% Now extract the stepnum from : %s.%10.10d.data
116			ic = findstr(filin,fil)+length(fil)+1; i = 1; clear stepnum
117			while filin(ic) ~= '.'
118			stepnum(i) = filin(ic); i = i + 1; ic = ic + 1;
119			end
120			ID = str2num(stepnum);
121			TIME(il,:) = datestr(datenum(1992,1,1)+ID*dt/60/60/24,'yyyymmddHHMM');
122			end
123			nt = size(TIME,1);
124
125			% %%%%%%%%%%%%%%%
126			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over time
127			for it = 1 : nt
128			snapshot = TIME(it,:);
129			ID = 606024/dt*( datenum(snapshot,'yyyymmddHHMM') - datenum(1992,1,1) );
130			filin = ll(it).name;
131			disp('')
132			disp(strcat('Processing: ',fil,'//',snapshot))
133			dirout = strcat(patho,sla,TIME(it,:),sla);
134			filout = sprintf('%s.%s.nc',otab{ifield,1},sub_name);
135
136			if ~exist(strcat(dirout,sla,filout),'file') % File already exists ?
137
138			%%%% READ THE FILE
139			switch otab{ifield,6}
140			case 4, flt = 'float32';
141			case 8, flt = 'float64';
142			end
143			t0 = clock;
144			if findstr(filin,'.gz') % Gunzip file, special care !
145			disp('\|----> Find a file with gz extension, work on uncompressed file ...')
146
147			% 1: copy the filename with it path into gunzip_1_file.txt
148			fid1 = fopen('gunzip_1_file.txt','w');
149			fprintf(fid1,'%s',strcat(pathi,sla,fild,sla,filin));fclose(fid1);
150
151			% 2: uncompress the file into a temporary folder:
152			disp('\|------> uncompressing the file ...')
153			! ./gunzip_1_file.bat
154			disp(strcat('\|--------: ',num2str(etime(t0,clock))))
155
156			% 3: Read the uncompress file:
157			disp('\|--> reading it ...')
158			C = readrec_cs510(strcat('gunzip_1_file',sla,'tempo.data'),LIMITS(6),flt);
159			disp(strcat('\|----: ',num2str(etime(t0,clock))))
160
161			% 4: Suppress it
162			! \rm ./gunzip_1_file/tempo.data
163
164			else % Simply read the file:
165			disp('\|--> reading it ...')
166			C = readrec_cs510(strcat(pathi,sla,fild,sla,filin),LIMITS(6),flt);
167			disp(strcat('\|----: ',num2str(etime(t0,clock))))
168			end
169
170			%%%% RESTRICT TO SUBDOMAIN
171			disp('\|--> get subdomain ...')
172			% Restrict vertical to subdomain:
173			if LIMITS(5) ~= 1
174			disp('\|----> vertical ...');
175			C = C(:,:,LIMITS(5):end);
176			end
177			% Clean the field:
178			C(find(C==0)) = NaN;
179			% Move the field into lat/lon grid:
180			disp('\|----> Move to lat/lon grid ...');
181			C = cube2latlon_fast(del,C);
182			% And then restrict horizontal to subdomain:
183			disp('\|----> horizontal ...');
184			C = C(max(find(XI<=LIMITS(1))):max(find(XI<=LIMITS(2))),...
185			max(find(YI<=LIMITS(3))):max(find(YI<=LIMITS(4))),:);
186
187
188			%%%% RECORD
189			disp('\|--> record netcdf file ...')
190			fid1 = fopen('inprogress.txt','w');
191			fprintf(fid1,'%s',strcat(dirout,sla,filout));fclose(fid1);
192
193			if 1 % Realy want to record ?
194
195			if ~exist(dirout,'dir')
196			mkdir(dirout);
197			end
198
199
200			nc = netcdf('inprogress.nc','clobber');
201
202			nc('X') = length(xi);
203			nc('Y') = length(yi);
204			nc('Z') = length(zi);
205
206			nc{'X'}='X';
207			nc{'Y'}='Y';
208			nc{'Z'}='Z';
209
210			nc{'X'}.uniquename='X';
211			nc{'X'}.long_name='longitude';
212			nc{'X'}.gridtype=ncint(0);
213			nc{'X'}.units='degrees_east';
214			nc{'X'}(:) = xi;
215
216			nc{'Y'}.uniquename='Y';
217			nc{'Y'}.long_name='latitude';
218			nc{'Y'}.gridtype=ncint(0);
219			nc{'Y'}.units='degrees_north';
220			nc{'Y'}(:) = yi;
221
222			nc{'Z'}.uniquename='Z';
223			nc{'Z'}.long_name='depth';
224			nc{'Z'}.gridtype=ncint(0);
225			nc{'Z'}.units='m';
226			nc{'Z'}(:) = zi;
227
228			ncid = fil;
229			nc{ncid}={'Z' 'Y' 'X'};
230			nc{ncid}.missing_value = ncdouble(NaN);
231			nc{ncid}.FillValue_ = ncdouble(0.0);
232			nc{ncid}(:,:,:) = permute(C,[3 2 1]);
233			nc{ncid}.units = otab{ifield,5};
234
235			close(nc);
236			! ./inprogress.bat
237
238			end %if 1/0 want to record ?
239			disp(strcat('\|--: ',num2str(etime(t0,clock))))
240
241			else
242			disp(strcat('\|--> Skip file (already done):',dirout,sla,filout))
243			end %if %file exist
244
245			end %for it
246			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END Loop over time
247			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
248
249			end %if it
250			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END Loop over variables to read
251			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%