/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_IO/process2nctiles.m
ViewVC logotype

Diff of /MITgcm_contrib/gael/matlab_class/gcmfaces_IO/process2nctiles.m

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

revision 1.3 by gforget, Sun Feb 2 01:47:38 2014 UTC revision 1.13 by gforget, Wed Mar 16 15:02:19 2016 UTC
# Line 1  Line 1 
1  function []=process2nctiles(dirModel,fileModel,fldModel);  function []=process2nctiles(dirModel,fileModel,fldModel,tileSize);
2  %process2nctiles(dirModel);  %process2nctiles(dirModel);
3  %object : convert MITgcm binary output to netcdf files (tiled)  %object : convert MITgcm binary output to netcdf files (tiled)
4  %inputs : dirModel is the MITgcm run directory  %inputs : dirModel is the MITgcm run directory
# Line 9  function []=process2nctiles(dirModel,fil Line 9  function []=process2nctiles(dirModel,fil
9  %           By default : all variables in e.g. 'state_2d_set1*'  %           By default : all variables in e.g. 'state_2d_set1*'
10  %           files will be processed, and writen individually to  %           files will be processed, and writen individually to
11  %           nctiles (tiled netcdf) that will be located in 'nctiles/'  %           nctiles (tiled netcdf) that will be located in 'nctiles/'
12  %         fldModel (optional) can be specifiec as e.g. 'ETAN'  %         fldModel (by default []) can be specified (as e.g. 'ETAN')
13    %            when fldModel is empty, all fields are processed
14    %         tileSize (optional) is e.g. [90 90] (by default tiles=faces)
15  %output : (netcdf files)  %output : (netcdf files)
16    
17  gcmfaces_global;  gcmfaces_global;
18    
19    %listFiles={'state_2d_set1','state_2d_set2','state_3d_set1','state_3d_set2'};
20    %listFiles={'trsp_3d_set1','trsp_3d_set2','trsp_3d_set3'};
21    %for ff=1:length(listFiles); process2nctiles('iter12/',listFiles{ff},[],[90 90]); end;
22    
23    %replace time series with monthly climatology?
24    doClim=0;
25    
26  %directory names  %directory names
27  listDirs={'STATE/','TRSP/'};  listDirs={'STATE/','TRSP/'};%BUDG?
28  filAvailDiag=[dirModel 'available_diagnostics.log'];  filAvailDiag=[dirModel 'available_diagnostics.log'];
29  filReadme=[dirModel 'README'];  filReadme=[dirModel 'README'];
30  dirOut=[dirModel 'nctiles/'];  dirOut=[dirModel 'nctiles_tmp/'];
31    %dirOut=[dirModel 'nctiles_post_tmp/'];
32  if ~isdir(dirOut); mkdir(dirOut); end;  if ~isdir(dirOut); mkdir(dirOut); end;
33    
34  %seacrh in subdirectories  %search in subdirectories
35  subDir=[];  subDir=[];
36    diagsDir='diags/';
37    %diagsDir='diags_post/';
38    %diagsDir='diags_interp/';
39  for ff=1:length(listDirs);  for ff=1:length(listDirs);
40  tmp1=dir([dirModel 'diags/' listDirs{ff} fileModel '*']);  tmp1=dir([dirModel diagsDir listDirs{ff} fileModel '*']);
41  if ~isempty(tmp1); subDir=listDirs{ff}; end;  if ~isempty(tmp1); subDir=listDirs{ff}; end;
42  end;  end;
43    
44    if isempty(subDir);
45    tmp1=dir([dirModel diagsDir fileModel '/' fileModel '*']);
46    if ~isempty(tmp1); subDir=[fileModel '/']; end;
47    end;
48    
49  if isempty(subDir);  if isempty(subDir);
50    error(['file ' fileModel ' was not found']);    error(['file ' fileModel ' was not found']);
51  else;  else;
52    dirIn=[dirModel 'diags/' subDir];    dirIn=[dirModel diagsDir subDir];
53    nn=length(dir([dirIn fileModel '*data']));    nn=length(dir([dirIn fileModel '*data']));
54    fprintf('%s (%d files) was found in \n %s \n',fileModel,nn,dirIn);    fprintf('%s (%d files) was found in \n %s \n',fileModel,nn,dirIn);
55  end;  end;
56    
57  %set list of variables to process  %set list of variables to process
58  if ~isempty(who('fldModel'));  if ~isempty(fldModel);
59     if ischar(fldModel); listFlds={fldModel};     if ischar(fldModel); listFlds={fldModel};
60     else; listFlds=fldModel;     else; listFlds=fldModel;
61     end;     end;
# Line 46  else; Line 64  else;
64     listFlds=meta.fldList;     listFlds=meta.fldList;
65  end;  end;
66    
67    %determine map of tile indices (by default tiles=faces)
68    if isempty(whos('tileSize'));
69      tileNo=mygrid.XC;
70      for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
71    else;
72      tileNo=gcmfaces_loc_tile(tileSize(1),tileSize(2));
73    end;
74    
75  %now do the actual processing  %now do the actual processing
76  for vv=1:length(listFlds);  for vv=1:length(listFlds);
77  nameDiag=deblank(listFlds{vv})  nameDiag=deblank(listFlds{vv})
# Line 58  if length(irec)~=1; error('field not in Line 84  if length(irec)~=1; error('field not in
84  %read time series  %read time series
85  myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);  myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);
86    
87  %ancilliary fields for netcdf file  %replace time series with monthly climatology
88    if doClim; myDiag=compClim(myDiag); end;
89    
90    %set ancilliary time variable
91  nn=length(size(myDiag{1}));  nn=length(size(myDiag{1}));
92  tim=[1:size(myDiag{1},nn)];  nn=size(myDiag{1},nn);
93  descr=nameDiag;  %tim=[1:nn];
94    tim=[1992*ones(nn,1) [1:nn]' 15*ones(nn,1)];
95    tim=datenum(tim)-datenum([1992 1 0]);
96    timUnits='days since 1992-1-1 0:0:0';
97    
98    %get time step axis
99    [listTimes]=diags_list_times({dirIn},{fileModel});
100    
101  %get units and long name from available_diagnostics.log  %get units and long name from available_diagnostics.log
102  [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);  [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
# Line 70  descr=nameDiag; Line 105  descr=nameDiag;
105  [rdm]=read_readme(filReadme);  [rdm]=read_readme(filReadme);
106  disp(rdm');  disp(rdm');
107    
108  %convert to MITgcm format (90x1170 array)  %set output directory/file name
109  myFile=[dirOut nameDiag];  myFile=[dirOut nameDiag];%first instance is for subdirectory name
110    if ~isdir(myFile); mkdir(myFile); end;
111    myFile=[myFile filesep nameDiag];%second instance is for file name base
112    
113  %get grid params  %get grid params
114  [grid_diag]=set_grid_diag(avail_diag);  [grid_diag]=set_grid_diag(avail_diag);
115    
116  %create netcdf file using write2nctiles (works only with old matlab, thus far ...)  %apply mask, and convert to land mask
117    if ~isempty(mygrid.RAC);
118      msk=grid_diag.msk;
119      if length(size(myDiag{1}))==3;
120        msk=repmat(msk(:,:,1),[1 1 size(myDiag{1},3)]);
121      else;
122        msk=repmat(msk,[1 1 1 size(myDiag{1},4)]);
123      end;
124      myDiag=myDiag.*msk;
125      clear msk;
126      %
127      land=isnan(grid_diag.msk);
128    end;
129    
130    %set 'coord' attribute
131    if avail_diag.nr~=1;
132      coord='lon lat dep tim';
133    else;
134      coord='lon lat tim';
135    end;
136    
137    %replace time series with monthly climatology
138    if doClim;
139      listTimes=listTimes(1:12);
140      timUnits='days since year-1-1 0:0:0';
141      avail_diag.longNameDiag=[avail_diag.longNameDiag ' (climatology) '];
142    end;
143    
144    %create netcdf file using write2nctiles
145  doCreate=1;  doCreate=1;
146  dimOut=write2nctiles(myFile,myDiag,doCreate,...  dimlist=write2nctiles(myFile,myDiag,doCreate,{'tileNo',tileNo},...
147      {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...      {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
148      {'units',avail_diag.units},{'descr',descr},{'rdm',rdm});      {'units',avail_diag.units},{'descr',nameDiag},{'coord',coord},{'rdm',rdm});
149    
150  %prepare to add fields  %determine relevant dimensions
151  doCreate=0;  for ff=1:length(dimlist);
152      dim.tim{ff}={dimlist{ff}{1}};
153  %determine relevant dimensions (note the reverse order)    dim.twoD{ff}={dimlist{ff}{end-1:end}};
154  for ff=1:mygrid.nFaces;    if avail_diag.nr~=1;
155    dimtim{ff}={dimOut{ff}{1}};      dim.threeD{ff}={dimlist{ff}{end-2:end}};
156    dim2d{ff}={dimOut{ff}{end-1:end}};      dim.dep{ff}={dimlist{ff}{end-2}};
   if avail_diag.nr~=1;  
     dimmsk{ff}={dimOut{ff}{end-2:end}};  
     dimdep{ff}={dimOut{ff}{end-2}};  
157    else;    else;
158      dimmsk{ff}=dim2d{ff};      dim.threeD{ff}=dim.twoD{ff};
159      dimdep{ff}=[];      dim.dep{ff}=[];
160    end;    end;
161  end;  end;
162    
163    %prepare to add fields
164    doCreate=0;
165    
166  %now add fields  %now add fields
167  write2nctiles(myFile,tim,doCreate,{'fldName','month'},...  write2nctiles(myFile,grid_diag.lon,doCreate,{'tileNo',tileNo},...
168    {'longName','month index starting Jan. 1992'},...    {'fldName','lon'},{'units','degrees_east'},{'dimIn',dim.twoD});
169    {'units',''},{'dimIn',dimtim});  write2nctiles(myFile,grid_diag.lat,doCreate,{'tileNo',tileNo},...
170  write2nctiles(myFile,grid_diag.lon,doCreate,...    {'fldName','lat'},{'units','degrees_north'},{'dimIn',dim.twoD});
   {'fldName','longitude'},{'units',''},{'dimIn',dim2d});  
 write2nctiles(myFile,grid_diag.lat,doCreate,...  
   {'fldName','latitude'},{'units',''},{'dimIn',dim2d});  
 write2nctiles(myFile,grid_diag.msk,doCreate,...  
   {'fldName','landmask'},{'units',''},...  
   {'longName','land mask'},{'dimIn',dimmsk});  
 write2nctiles(myFile,grid_diag.RAC,doCreate,...  
   {'fldName','cellarea'},{'units','m^2'},...  
   {'longName','grid cell area'},{'dimIn',dim2d});  
171  if isfield(grid_diag,'dep');  if isfield(grid_diag,'dep');
172      write2nctiles(myFile,grid_diag.dep,doCreate,...      write2nctiles(myFile,grid_diag.dep,doCreate,{'tileNo',tileNo},...
173        {'fldName','depth'},{'units','m'},{'dimIn',dimdep});        {'fldName','dep'},{'units','m'},{'dimIn',dim.dep});
174      write2nctiles(myFile,grid_diag.dz,doCreate,...  end;
175        {'fldName','cellthick'},{'units','m'},{'dimIn',dimdep});  write2nctiles(myFile,tim,doCreate,{'tileNo',tileNo},...
176      {'fldName','tim'},{'longName','time'},...
177      {'units',timUnits},{'dimIn',dim.tim});
178    if ~isempty(mygrid.RAC);
179      write2nctiles(myFile,listTimes,doCreate,{'tileNo',tileNo},...
180        {'fldName','timstep'},{'longName','final time step number'},...
181        {'units','1'},{'dimIn',dim.tim});
182      write2nctiles(myFile,grid_diag.msk,doCreate,{'tileNo',tileNo},...
183        {'fldName','land'},{'units','1'},{'longName','land mask'},{'dimIn',dim.threeD});
184      write2nctiles(myFile,grid_diag.RAC,doCreate,{'tileNo',tileNo},...
185        {'fldName','area'},{'units','m^2'},{'longName','grid cell area'},{'dimIn',dim.twoD});
186      if isfield(grid_diag,'dz');
187        write2nctiles(myFile,grid_diag.dz,doCreate,{'tileNo',tileNo},...
188          {'fldName','thic'},{'units','m'},{'dimIn',dim.dep});
189      end;
190  end;  end;
191    
192    clear myDiag;
193    
194  end;%for vv=1:length(listFlds);  end;%for vv=1:length(listFlds);
195    
196  function [meta]=read_meta(fileName);  function [meta]=read_meta(fileName);
197    
198  %read meta file  %read meta file
199  tmp1=dir([fileName '*.meta']); tmp1=tmp1(1).name;  tmp1=dir([fileName '*.meta']); tmp1=tmp1(1).name;
200  tmp2=strfind(fileName,'/');  tmp2=strfind(fileName,filesep);
201  if ~isempty(tmp2); tmp2=tmp2(end); else; tmp2=0; end;  if ~isempty(tmp2); tmp2=tmp2(end); else; tmp2=0; end;
202  tmp1=[fileName(1:tmp2) tmp1]; fid=fopen(tmp1);  tmp1=[fileName(1:tmp2) tmp1]; fid=fopen(tmp1);
203  while 1;  while 1;
# Line 236  if avail_diag.nr~=1; Line 308  if avail_diag.nr~=1;
308      grid_diag.dep=reshape(grid_diag.dep,[1 1 avail_diag.nr]);      grid_diag.dep=reshape(grid_diag.dep,[1 1 avail_diag.nr]);
309      grid_diag.dz=reshape(grid_diag.dz,[1 1 avail_diag.nr]);      grid_diag.dz=reshape(grid_diag.dz,[1 1 avail_diag.nr]);
310  end;  end;
311    
312    %%replace time series with monthly climatology
313    function [FLD]=compClim(fld);
314    
315    gcmfaces_global;
316    
317    ndim=length(size(fld{1}));
318    nyear=size(fld{1},ndim)/12;
319    
320    if ndim==3; FLD=NaN*fld(:,:,1:12); end;
321    if ndim==4; FLD=NaN*fld(:,:,:,1:12); end;
322    
323    for mm=1:12;
324    if ndim==3; FLD(:,:,mm)=mean(fld(:,:,mm:12:12*nyear),ndim); end;
325    if ndim==4; FLD(:,:,:,mm)=mean(fld(:,:,:,mm:12:12*nyear),ndim); end;
326    end;
327    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22