Nctilesdev (#3)

* - add dirDiags as input argument (interp2nctiles.m, process2interp.m)
- add dimlist, dimname, and clmbnds optional arguments (write2nctiles.m)
- revise handling of directories and time loop (process2interp.m, process2nctiles.m)
- rename dirModel,fileModel,fldModel as dirDiags,fileDiags,selectFld (process2nctiles.m)
- revise handling of doClim option and time axis using clmbnds (process2nctiles.m, write2nctiles.m)
- add possibility to rename variables on the fly using filRename  (process2nctiles.m)
- rename coordinates as 't','k','j','i' (previously 'i1','i2','i3','i4') via dimlist,
  dimname being passed as arguments to write2nctiles.m (process2nctiles.m).
- set 'xtype' to 'float' for main variable to reduce file sizes (process2nctiles.m).
- revise handling of dimension, file, and variable names (write2nctiles.m).

* Changes in grid_load.m (XW etc.), gcmfaces_loc_tile.m (tile numbering), nctiles I/O (dimensions etc.).
- grid_load.m: (re)assign RF(nr+1) and DRC(nr+1) if needed; compute XW/YW, XS/YS from XC/YC.
- gcmfaces_loc_tile.m: switch to ordering convention that is consistent with MITgcm/pkg/exch2.
- read_nctiles.m: add nctiles_old_tile_order switch so that the user can revert to the old ordering if needed (only matters for smaller tile sizes).
- gcmfaces_IO/write2nctiles.m: add itile, ntile attributes (their presence will notably indicate that the new tile ordering convention was used); revisit the handling of dimlist in the doCreate 1/0 phases (needed after previous commit).
- process2nctiles.m: fix bug introduced in previous commit (tim in doClim case & handling of filRename); revisit set_grid_diag function to use predefined dimension names that distinguish C,W,S and C,L,U points (e.g., {'t','kc','jc','ic'} for tracer fields) and set the longitude/latitude variables accordingly (e.g., use XW/YW for velocity components).
- interp2nctiles.m: comment out adhoc DRC, DRF, mskC definitions to avoid putting these in the netcdf files.
- write2nctiles.m: add itile, mtile attributes; clean up the two stage handling of dimlist, dimname, dimvec.

* - interp2nctiles.m : update to match revised process2nctiles.m argument handling.
- process2interp.m : go directly from mds output to interpolation (rather than via
  nctiles first); revise handling of input arguments (see help section for detail).
- process2nctiles.m : revise handling of input arguments (see help section for
  detail); use rdmds_meta.m in place of embeded duplicate (read_meta function).
- rdmds_search_subdirs.m (new) : searches for fileDiags within the subdirectories
  of dirDiags (now used in process2interp.m and process2nctiles.m)

* - process2UVSTAR.m (new) computes bolus velocity from GM_PsiX,Y
- process2UEVN.m get dirDiags,fileDiags via function arguments (as now done in process*.m),
  use rdmds_search_subdirs (to identify subDir), and calc_UV_zonmer (no masking unlike before).

* changes related to rdmds_meta.m and process2interp.m
- deal with adding '*' to file name within rdmds_meta.m (not in caller)
- add possibility to rely on interp_precomputed.mat in process2interp.m

* improve treatment of missing_value, _FillValue, scale_factor (new), and add_offset (new)

Author: gaelforget

gcmfaces_IO/grid_load.m

@@ -149,6 +149,11 @@
     end;
 end;
 
+%if needed then (re)set RF(nr+1) and DRC(nr+1):
+nr=length(mygrid.RC);
+if isnan(mygrid.RF(nr+1)); mygrid.RF(nr+1)=mygrid.RC(nr)-1/2*mygrid.DRF(nr); end;
+if length(mygrid.DRC)==nr||isnan(mygrid.DRC(nr+1)); mygrid.DRC(nr+1)=mygrid.DRC(nr)/2; end;
+
 %grid orientation
 if mygrid.memoryLimit<2;
     list0={'AngleCS','AngleSN'};
@@ -198,13 +203,50 @@
 %grid masks
 if mygrid.memoryLimit<1;
     mygrid.hFacCsurf=mygrid.hFacC;
-    for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
+    for ff=1:mygrid.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
 
     mskC=mygrid.hFacC; mskC(mskC==0)=NaN; mskC(mskC>0)=1; mygrid.mskC=mskC;
     mskW=mygrid.hFacW; mskW(mskW==0)=NaN; mskW(mskW>0)=1; mygrid.mskW=mskW;
     mskS=mygrid.hFacS; mskS(mskS==0)=NaN; mskS(mskS>0)=1; mygrid.mskS=mskS;
 end;
 
+%missing velocity point locations
+if mygrid.memoryLimit<2;
+    %1) average C point locations
+    tmpXC=exch_T_N(mygrid.XC); tmpYC=exch_T_N(mygrid.YC);
+    mygrid.XW=NaN*mygrid.XC; mygrid.YW=NaN*mygrid.YC;
+    mygrid.XS=NaN*mygrid.XC; mygrid.YS=NaN*mygrid.YC;
+    for ff=1:mygrid.nFaces;
+       tmp1=tmpXC{ff}(1:end-2,2:end-1);
+       tmp2=tmpXC{ff}(2:end-1,2:end-1); 
+       tmp2(tmp2-tmp1>180)=tmp2(tmp2-tmp1>180)-360;
+       tmp2(tmp1-tmp2>180)=tmp2(tmp1-tmp2>180)+360;
+       mygrid.XW{ff}=(tmp1+tmp2)/2;
+       %
+       tmp1=tmpXC{ff}(2:end-1,1:end-2);
+       tmp2=tmpXC{ff}(2:end-1,2:end-1);
+       tmp2(tmp2-tmp1>180)=tmp2(tmp2-tmp1>180)-360;
+       tmp2(tmp1-tmp2>180)=tmp2(tmp1-tmp2>180)+360;
+       mygrid.XS{ff}=(tmp1+tmp2)/2;
+       %
+       tmp1=tmpYC{ff}(1:end-2,2:end-1);
+       tmp2=tmpYC{ff}(2:end-1,2:end-1);
+       mygrid.YW{ff}=(tmp1+tmp2)/2;
+       %
+       tmp1=tmpYC{ff}(2:end-1,1:end-2);
+       tmp2=tmpYC{ff}(2:end-1,2:end-1);
+       mygrid.YS{ff}=(tmp1+tmp2)/2;
+    end;
+    %2) fix lonitude range:
+    tmp1=convert2gcmfaces(mygrid.XC);
+    test1=~isempty(find(tmp1(:)<0));
+    if test1; Xmax=180; Xmin=-180; else; Xmax=360; Xmin=0; end;
+    mygrid.XS(mygrid.XS<Xmin)=mygrid.XS(mygrid.XS<Xmin)+360;
+    mygrid.XS(mygrid.XS>Xmax)=mygrid.XS(mygrid.XS>Xmax)-360;
+    mygrid.XW(mygrid.XW<Xmin)=mygrid.XW(mygrid.XW<Xmin)+360;
+    mygrid.XW(mygrid.XW>Xmax)=mygrid.XW(mygrid.XW>Xmax)-360;
+end;
+
 %zonal mean and sections needed for transport computations
 % if mygrid.memoryLimit<1;
 %     if ~isfield(mygrid,'mygrid.LATS_MASKS');

gcmfaces_IO/interp2nctiles.m

@@ -1,5 +1,6 @@
-function []=interp2nctiles(indFiles);
-% INTERP2CNTILES creates netcdf files from interpolated 
+function []=interp2nctiles(dirDiags,listDo);
+% INTERP2CNTILES(dirDiags,listDo); 
+%   creates netcdf files from interpolated 
 %   fields that were created by process2interp. The 
 %   input vector process2interp specifies the file subset
 %   to be processed (1:length(listInterp) by default).
@@ -9,11 +10,7 @@
 
 gcmfaces_global; global mygrid_orig;
 
-dirModel='./';
-if isempty(mygrid_orig);
-grid_load; mygrid_orig=mygrid;
-cd(dirModel);
-end;
+if isempty(mygrid_orig); mygrid_orig=mygrid; end;
 
 lon=[-179.75:0.5:179.75]; lat=[-89.75:0.5:89.75];
 [lat,lon] = meshgrid(lat,lon);
@@ -26,22 +23,16 @@
 mygrid_latlon.YC=gcmfaces({lat});
 mygrid_latlon.RC=mygrid.RC;
 mygrid_latlon.RF=mygrid.RF;
-mygrid_latlon.DRC=mygrid.DRC;
-mygrid_latlon.DRF=mygrid.DRF;
-mygrid_latlon.mskC=1+0*repmat(mygrid_latlon.XC,[1 1  length(mygrid.RC)]);
-mygrid_latlon.RAC=[];
+%mygrid_latlon.DRC=mygrid.DRC;
+%mygrid_latlon.DRF=mygrid.DRF;
+%mygrid_latlon.mskC=1+0*repmat(mygrid_latlon.XC,[1 1  length(mygrid.RC)]);
 mygrid_latlon.gcm2facesFast=0;
 mygrid_latlon.facesExpand=[];
 
-[listInterp,listNot]=process2interp;
-if isempty(who('indFiles')); 
-  indFiles=[1:length(listInterp)];
-end;
-
 mygrid=mygrid_latlon;
-for ii=indFiles;
-tic; process2nctiles(dirModel,listInterp{ii},[]);
-fprintf(['DONE: ' listInterp{ii} ' (in ' num2str(toc) 's)\n']);
+for ii=1:length(listDo);
+tic; process2nctiles(dirDiags,listDo{ii},listDo{ii});
+fprintf(['DONE: ' listDo{ii} ' (in ' num2str(toc) 's)\n']);
 end;
 mygrid=mygrid_orig;
 

gcmfaces_IO/ncload.m

@@ -25,9 +25,17 @@
             bb=length(size(aa)); aa=permute(aa,[bb:-1:1]);
             %replace missing value with NaN
             [atts]=ncatts(f,varid);
-            if strcmp(atts,'missing_value')&isreal(aa);
+            if sum(strcmp(atts,'scale_factor'));
+                scale_factor=double(netcdf.getAtt(f,varid,'scale_factor'));
+                aa=scale_factor*double(aa);
+            end;            
+            if sum(strcmp(atts,'add_offset'));
+                add_offset=double(netcdf.getAtt(f,varid,'add_offset'));
+                aa=aa+add_offset;
+            end;            
+            if sum(strcmp(atts,'missing_value'))&isreal(aa);
                 spval = double(netcdf.getAtt(f,varid,'missing_value'));
-            elseif strcmp(atts,'_FillValue')&isreal(aa);
+            elseif sum(strcmp(atts,'_FillValue'))&isreal(aa);
                 spval = double(netcdf.getAtt(f,varid,'_FillValue'));
             else;
                 spval=[];

gcmfaces_IO/process2UEVN.m

@@ -1,64 +1,67 @@
-function []=process2UEVN(ff);
-% PROCESS2UEVN(ff) computes eastward/northward vector components
-%   and output result to 'diags_post_tmp/'
+function []=process2UEVN(dirDiags,fileDiags);
+% PROCESS2UEVN(dirDiags,fileDiags) computes eastward/northward vector components
+%   for vector fields in fileDiags and output result to 'diags_UEVN_tmp/'
 
 gcmfaces_global;
 
-%ff=1;
-dir0=[pwd filesep];
-dirIn=[dir0 'diags/'];
-dirOut=[dir0 'diags_post_tmp/'];
-if ~isdir(dirOut); mkdir(dirOut); end;
+gcmfaces_global;
+dirOut=[dirDiags filesep 'diags_UEVN_tmp' filesep];
 
-if ff==1;
-   dirIn = [dirIn 'TRSP/'];
-   dirOut = [dirOut 'TRSP/'];
-   fil = 'trsp_3d_set1';
+if strcmp(fileDiags,'trsp_3d_set1');
    pairsIn={{'UVELMASS','VVELMASS'}};
    pairsOut={{'EVELMASS','NVELMASS'}};
-elseif ff==2;
-   dirIn = [dirIn 'TRSP/'];
-   dirOut = [dirOut 'TRSP/'];
-   fil = 'trsp_3d_set2';
+elseif strcmp(fileDiags,'trsp_3d_set2');
    pairsIn={{'DFxE_TH ','DFyE_TH '},{'DFxE_SLT','DFyE_SLT'}};
    pairsOut={{'DFeE_TH ','DFnE_TH '},{'DFeE_SLT','DFnE_SLT'}};
    pairsIn={pairsIn{:},{'ADVx_TH ','ADVy_TH '},{'ADVx_SLT','ADVy_SLT'}};
    pairsOut={pairsOut{:},{'ADVe_TH ','ADVn_TH '},{'ADVe_SLT','ADVn_SLT'}};
-elseif ff==3;
-   dirIn = [dirIn 'STATE/'];
-   dirOut = [dirOut 'STATE/'];
-   fil= 'state_2d_set1';
+elseif strcmp(fileDiags,'state_2d_set1');
    pairsIn={{'DFxEHEFF','DFyEHEFF'},{'DFxESNOW','DFyESNOW'}};
    pairsOut={{'DFeEHEFF','DFnEHEFF'},{'DFeESNOW','DFnESNOW'}};
    pairsIn={pairsIn{:},{'ADVxHEFF','ADVyHEFF'},{'ADVxSNOW','ADVySNOW'}};
    pairsOut={pairsOut{:},{'ADVeHEFF','ADVnHEFF'},{'ADVeSNOW','ADVnSNOW'}};
    pairsIn={pairsIn{:},{'oceTAUX ','oceTAUY '},{'SIuice  ','SIvice  '}};
    pairsOut={pairsOut{:},{'oceTAUE ','oceTAUN '},{'EVELice ','NVELice '}};
+elseif strcmp(fileDiags,'star_trsp_3d_set1');
+   pairsIn={{'UVELSTAR','VVELSTAR'}};
+   pairsOut={{'EVELSTAR','NVELSTAR'}};
+else;
+   error('unknown fileDiags');
 end;
 
-listIn=dir([dirIn fil '*meta']);
+%% ======== PART 1 =======
+
+%search for fileDiags in subdirectories
+[subDir]=rdmds_search_subdirs(dirDiags,fileDiags);
+%read meta file to get list of variables
+[meta]=rdmds_meta([dirDiags subDir fileDiags]);
+
+%% ======== PART 2 =======
+
+listIn=dir([dirDiags subDir fileDiags '*meta']);
 for tt=1:length(listIn);
   disp([tt length(listIn)]);
   fldsOut=gcmfaces;
   listOut={};
   for pp=1:length(pairsIn);
     pIn=pairsIn{pp}; pOut=pairsOut{pp};
     filIn=listIn(tt).name(1:end-5);
-    metaIn=rdmds_meta([dirIn filIn]);
+    metaIn=rdmds_meta([dirDiags subDir filIn]);
     i1=find(strcmp(metaIn.fldList,pIn{1}));
     i2=find(strcmp(metaIn.fldList,pIn{2}));
     %[i1 i2]
     %
-    UX=rdmds2gcmfaces([dirIn filIn],'rec',i1);
-    VY=rdmds2gcmfaces([dirIn filIn],'rec',i2);
-    [UE,VN]=calc_UEVNfromUXVY(UX,VY);
+    UX=rdmds2gcmfaces([dirDiags subDir filIn],'rec',i1);
+    VY=rdmds2gcmfaces([dirDiags subDir filIn],'rec',i2);
+    [UE,VN]=calc_UV_zonmer(UX,VY);
+%   [UE,VN]=calc_UEVNfromUXVY(UX,VY);
     %store binary
     fldsOut=cat(3,fldsOut,UE);
     fldsOut=cat(3,fldsOut,VN);
     listOut={listOut{:},pOut{:}};
   end;
   %output binary file
-  filOut=['post_' filIn];
+  filOut=['zonmer_' filIn];
   tmp1=convert2gcmfaces(fldsOut);
   tmp1(isnan(tmp1))=0;
   if ~isdir(dirOut); mkdir(dirOut); end;

gcmfaces_IO/process2UVSTAR.m

@@ -0,0 +1,53 @@
+function []=process2UVSTAR(dirDiags,fileDiags);
+% PROCESS2UVstar computes bolus velocity vector components (U,V,WVELSTAR) from
+%   GM_PsiX,Y in fileDiags (trsp_3d_set1) and output result to 'diags_post_tmp/'
+
+gcmfaces_global;
+
+dirOut=[dirDiags filesep 'diags_UVstar_tmp' filesep];
+if ~isdir(dirOut); mkdir(dirOut); end;
+
+pairsIn={{'GM_PsiX ' 'GM_PsiY '}};
+pairsOut={{'UVELSTAR','VVELSTAR','WVELSTAR'}};
+
+%% ======== PART 1 =======
+
+%search for fileDiags in subdirectories
+[subDir]=rdmds_search_subdirs(dirDiags,fileDiags);
+%read meta file to get list of variables
+[meta]=rdmds_meta([dirDiags subDir fileDiags]);
+
+listIn=dir([dirDiags subDir fileDiags '*meta']);
+for tt=1:length(listIn);
+  disp([tt length(listIn)]);
+  fldsOut=gcmfaces;
+  listOut={};
+  for pp=1:length(pairsIn);
+    pIn=pairsIn{pp}; pOut=pairsOut{pp};
+    filIn=listIn(tt).name(1:end-5);
+    metaIn=rdmds_meta([dirDiags subDir filIn]);
+    i1=find(strcmp(metaIn.fldList,pIn{1}));
+    i2=find(strcmp(metaIn.fldList,pIn{2}));
+    %[i1 i2]
+    %
+    GM_PsiX=rdmds2gcmfaces([dirDiags subDir filIn],'rec',i1);
+    GM_PsiY=rdmds2gcmfaces([dirDiags subDir filIn],'rec',i2);
+    [fldUbolus,fldVbolus,fldWbolus]=calc_bolus(GM_PsiX,GM_PsiY);
+    %store binary
+    fldsOut=cat(3,fldsOut,fldUbolus);
+    fldsOut=cat(3,fldsOut,fldVbolus);
+    fldsOut=cat(3,fldsOut,fldWbolus);
+    listOut={listOut{:},pOut{:}};
+  end;
+  %output binary file
+  filOut=['star_' filIn];
+  tmp1=convert2gcmfaces(fldsOut);
+  tmp1(isnan(tmp1))=0;
+  if ~isdir(dirOut); mkdir(dirOut); end;
+  write2file([dirOut filOut '.data'],tmp1,32);
+  %create meta file
+  tmp2=size(tmp1); tmp2(end)=tmp2(end)/3;
+  write2meta([dirOut filOut '.data'],tmp2,32,listOut);
+end;
+
+