/[MITgcm]/MITgcm_contrib/mlosch/interp_llc/mdsiocompact.m
ViewVC logotype

Diff of /MITgcm_contrib/mlosch/interp_llc/mdsiocompact.m

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

revision 1.1 by mlosch, Thu May 3 21:07:21 2007 UTC revision 1.2 by mlosch, Fri May 4 14:53:00 2007 UTC
# Line 1  Line 1 
1  function b=mdsiocompact(a);  function b=mdsiocompact(a,flag);
2    % hack to convert to mdsio compact format    % hack to convert to and from mdsio compact format
3        
4    [nx,ny,nz,nt]=size(a);    if nargin < 2
5    if nx/4~=round(nx/4);      flag = 1;
     error('not a llc field');  
6    end    end
7      if flag
8        % convert from "nice" to compact format
9        [nx,ny,nz,nt]=size(a);
10        if nx/4~=round(nx/4);
11          error('not a llc field');
12        end
13    
14        n4=nx/4;
15        m=(ny-n4)/n4;
16        m4=m*n4;
17        for kt = 1:nt
18          for kz = 1:nz
19            for k=1:4
20              sides{k} = a([1:n4]+(k-1)*n4,1:m4,kz,kt);
21            end
22            cap = a([1:n4]+1*n4,m4+1:end,kz,kt);
23            % reformat
24            btmp = [sides{1} ...
25                    sides{2} ...
26                    cap ...
27                    reshape(rot90(sides{3},1),[n4 m4]) ...
28                    reshape(rot90(sides{4},1),[n4 m4]) ...
29                    cap*0];
30            b(:,:,kz,kt)=btmp;
31          end
32        end
33      else
34        % convert from compact to "nice" format
35        [nx,m,nz,nt]=size(a);
36        if m/nx~=round(m/nx);
37          error('not a compact format llc field');
38        end
39    
40    n4=nx/4;      fac = 1;
41    m=(ny-n4)/n4;      ntiles = m/nx;
42    m4=m*n4;      if mod(ntiles,2); % odd number of tiles
43    for kt = 1:nt        tilesPerFace = (ntiles-1)/4;
44      for kz = 1:nz      else
45        for k=1:4        tilesPerFace = (ntiles-2)/4;
46          sides{k} = a([1:n4]+(k-1)*n4,1:m4,kz,kt);      end
47        n4=tilesPerFace*nx;
48        for kt = 1:nt
49          for kz = 1:nz
50            for k=1:4
51              jj = [1:n4]+(k-1)*n4;
52              if k>2; jj = jj + nx; end
53              sides{k} = a(:,jj,kz,kt);
54            end
55            cap = rot90(a(:,[1:nx]+2*n4,kz,kt),1);
56            % reformat
57            btmp = [sides{1}; ...
58                    sides{2}; ...
59                    rot90(reshape(sides{3}(:),[n4 nx]),-1); ...
60                    rot90(reshape(sides{4}(:),[n4 nx]),-1) ];
61            b(:,:,kz,kt)=[btmp ...
62                          [cap; rot90(cap,-1)*abs(fac); ...
63                           rot90(cap,-2)*abs(fac); ...
64                           rot90(cap,-3)*abs(fac)]];
65        end        end
       cap = a([1:n4]+1*n4,m4+1:end,kz,kt);  
       % reformat  
       btmp = [sides{1} ...  
               sides{2} ...  
               cap ...  
               reshape(rot90(sides{3},1),[n4 m4]) ...  
               reshape(rot90(sides{4},1),[n4 m4]) ...  
               cap*0];  
       b(:,:,kz,kt)=btmp;  
66      end      end
67    end    end
68        

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

  ViewVC Help
Powered by ViewVC 1.1.22