/[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.1 by gforget, Thu Jan 30 17:05:12 2014 UTC revision 1.7 by gforget, Tue Feb 4 22:13:56 2014 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  fprintf('to do list : embed README in global attributes\n');  %listFiles={'state_2d_set1','state_2d_set2','state_3d_set1','state_3d_set2'};
20  fprintf('to do list : create tim variable\n');  %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  %directory names  %directory names
24  listDirs={'STATE/','TRSP/'};  listDirs={'STATE/','TRSP/'};
25  filAvailDiag=[dirModel 'available_diagnostics.log'];  filAvailDiag=[dirModel 'available_diagnostics.log'];
26    filReadme=[dirModel 'README'];
27  dirOut=[dirModel 'nctiles/'];  dirOut=[dirModel 'nctiles/'];
28  if ~isdir(dirOut); mkdir(dirOut); end;  if ~isdir(dirOut); mkdir(dirOut); end;
29    
# Line 39  else; Line 43  else;
43  end;  end;
44    
45  %set list of variables to process  %set list of variables to process
46  if ~isempty(who('fldModel'));  if ~isempty(fldModel);
47     if ischar(fldModel); listFlds={fldModel};     if ischar(fldModel); listFlds={fldModel};
48     else; listFlds=fldModel;     else; listFlds=fldModel;
49     end;     end;
# Line 48  else; Line 52  else;
52     listFlds=meta.fldList;     listFlds=meta.fldList;
53  end;  end;
54    
55    %determine map of tile indices (by default tiles=faces)
56    if isempty(whos('tileSize'));
57      tileNo=mygrid.XC;
58      for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
59    else;
60      tileNo=gcmfaces_loc_tile(tileSize(1),tileSize(2));
61    end;
62    
63  %now do the actual processing  %now do the actual processing
64  for vv=1:length(listFlds);  for vv=1:length(listFlds);
65  nameDiag=deblank(listFlds{vv})  nameDiag=deblank(listFlds{vv})
# Line 60  if length(irec)~=1; error('field not in Line 72  if length(irec)~=1; error('field not in
72  %read time series  %read time series
73  myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);  myDiag=rdmds2gcmfaces([dirIn fileModel '*'],NaN,'rec',irec);
74    
75  %ancilliary fields for netcdf file  %set ancilliary time variable
76  tim=[1:size(myDiag{1},3)];  nn=length(size(myDiag{1}));
77  timName='month index';  nn=size(myDiag{1},nn);
78  descr=nameDiag;  %tim=[1:nn];
79    tim=[1992*ones(nn,1) [1:nn]' 15*ones(nn,1)];
80    tim=datenum(tim)-datenum([1992 1 0]);
81    timUnits='days since 1992-1-1 0:0:0';
82    
83    %get time step axis
84    [listTimes]=diags_list_times({[dirModel 'diags/STATE/']},{'state_2d_set1'});
85    
86  %get units and long name from available_diagnostics.log  %get units and long name from available_diagnostics.log
87  [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);  [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
88    
89    %get description of estimate from README
90    [rdm]=read_readme(filReadme);
91    disp(rdm');
92    
93  %convert to MITgcm format (90x1170 array)  %convert to MITgcm format (90x1170 array)
94  myFile=[dirOut nameDiag];  myFile=[dirOut nameDiag];%first instance is for subdirectory name
95    if ~isdir(myFile); mkdir(myFile); end;
96    myFile=[myFile '/' nameDiag];%second instance is for file name base
97    
98  %get grid params  %get grid params
99  [grid_diag]=set_grid_diag(avail_diag);  [grid_diag]=set_grid_diag(avail_diag);
100    
101    %apply mask, and convert to land mask
102    msk=grid_diag.msk;
103    if length(size(myDiag{1}))==3;
104      msk=repmat(msk,[1 1 size(myDiag{1},3)]);
105    end;
106    for kk=1:size(myDiag{1},4);
107      myDiag(:,:,:,kk)=myDiag(:,:,:,kk).*msk;
108    end;
109    clear msk;
110    %
111    land=isnan(grid_diag.msk);
112    
113    %set 'coord' attribute
114    if avail_diag.nr~=1;
115      coord='lon lat dep tim';
116    else;
117      coord='lon lat tim';
118    end;
119    
120  %create netcdf file using write2nctiles (works only with old matlab, thus far ...)  %create netcdf file using write2nctiles (works only with old matlab, thus far ...)
121  doCreate=1;  doCreate=1;
122  dimOut=write2nctiles(myFile,myDiag,doCreate,...  dimlist=write2nctiles(myFile,myDiag,doCreate,{'tileNo',tileNo},...
123      {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...      {'fldName',nameDiag},{'longName',avail_diag.longNameDiag},...
124      {'units',avail_diag.units},{'descr',descr});      {'units',avail_diag.units},{'descr',nameDiag},{'coord',coord},{'rdm',rdm});
   
 %prepare to add fields  
 doCreate=0;  
125    
126  %determine relevant dimensions (note the reverse order)  %determine relevant dimensions
127  for ff=1:mygrid.nFaces;  for ff=1:length(dimlist);
128    dim2d{ff}={dimOut{ff}{end-1:end}};    dim.tim{ff}={dimlist{ff}{1}};
129    if avail_diag.nr~=1;    dim.twoD{ff}={dimlist{ff}{end-1:end}};
130      dimmsk{ff}={dimOut{ff}{end-2:end}};    if avail_diag.nr~=1;
131      dimdep{ff}={dimOut{ff}{end-2}};      dim.threeD{ff}={dimlist{ff}{end-2:end}};
132        dim.dep{ff}={dimlist{ff}{end-2}};
133    else;    else;
134      dimmsk{ff}=dim2d{ff};      dim.threeD{ff}=dim.twoD{ff};
135      dimdep{ff}=[];      dim.dep{ff}=[];
136    end;    end;
137  end;  end;
138    
139    %prepare to add fields
140    doCreate=0;
141    
142  %now add fields  %now add fields
143  write2nctiles(myFile,grid_diag.lon,doCreate,...  write2nctiles(myFile,tim,doCreate,{'tileNo',tileNo},...
144    {'fldName','longitude'},{'units',''},{'dimIn',dim2d});    {'fldName','tim'},{'longName','time'},...
145  write2nctiles(myFile,grid_diag.lat,doCreate,...    {'units',timUnits},{'dimIn',dim.tim});
146    {'fldName','latitude'},{'units',''},{'dimIn',dim2d});  write2nctiles(myFile,listTimes,doCreate,{'tileNo',tileNo},...
147  write2nctiles(myFile,grid_diag.msk,doCreate,...    {'fldName','step'},{'longName','final time step number'},...
148    {'fldName','landmask'},{'units',''},...    {'units','1'},{'dimIn',dim.tim});
149    {'longName','land mask'},{'dimIn',dimmsk});  write2nctiles(myFile,grid_diag.lon,doCreate,{'tileNo',tileNo},...
150  write2nctiles(myFile,grid_diag.RAC,doCreate,...    {'fldName','lon'},{'units','degrees_east'},{'dimIn',dim.twoD});
151    {'fldName','cellarea'},{'units','m^2'},...  write2nctiles(myFile,grid_diag.lat,doCreate,{'tileNo',tileNo},...
152    {'longName','grid cell area'},{'dimIn',dim2d});    {'fldName','lat'},{'units','degrees_north'},{'dimIn',dim.twoD});
153    write2nctiles(myFile,grid_diag.msk,doCreate,{'tileNo',tileNo},...
154      {'fldName','land'},{'units','1'},{'longName','land mask'},{'dimIn',dim.threeD});
155    write2nctiles(myFile,grid_diag.RAC,doCreate,{'tileNo',tileNo},...
156      {'fldName','area'},{'units','m^2'},{'longName','grid cell area'},{'dimIn',dim.twoD});
157  if isfield(grid_diag,'dep');  if isfield(grid_diag,'dep');
158      write2nctiles(myFile,grid_diag.dep,doCreate,...      write2nctiles(myFile,grid_diag.dep,doCreate,{'tileNo',tileNo},...
159        {'fldName','depth'},{'units','m'},{'dimIn',dimdep});        {'fldName','dep'},{'units','m'},{'dimIn',dim.dep});
160      write2nctiles(myFile,grid_diag.dz,doCreate,...      write2nctiles(myFile,grid_diag.dz,doCreate,{'tileNo',tileNo},...
161        {'fldName','cellthick'},{'units','m'},{'dimIn',dimdep});        {'fldName','thic'},{'units','m'},{'dimIn',dim.dep});
162  end;  end;
163    
164  end;%for vv=1:length(listFlds);  end;%for vv=1:length(listFlds);
# Line 144  if ~isempty(who('timeStepNumber'));  met Line 192  if ~isempty(who('timeStepNumber'));  met
192    
193  %%  %%
194    
195    function [rdm]=read_readme(filReadme);
196    
197    gcmfaces_global;
198    
199    rdm=[];
200    
201    fid=fopen(filReadme,'rt');
202    while ~feof(fid);
203        nn=length(rdm);
204        rdm{nn+1} = fgetl(fid);
205    end;
206    fclose(fid);
207    
208    %%
209    
210  function [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);  function [avail_diag]=read_avail_diag(filAvailDiag,nameDiag);
211    
212  gcmfaces_global;  gcmfaces_global;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22