Warning, /utils/matlab/cs_grid/rotate_uv2uvEN.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 6b4f212f on 2024-06-05 16:55:17 UTC
4582449b53 Jean*0001 function [uE,vN,msk] = rotate_uv2uvEN(u,v,AngleCS,AngleSN,Grid,maskW,maskS)
a8e064706f Jean*0002 % [uE,vN,msk] = rotate_uv2uvEN(u,v,AngleCS,AngleSN,[Grid,maskW,maskS])
5cea1a7045 Dani*0003 %
0004 % Rotate cube sphere U and V vector components to east-west (uE) and
0005 % north-south (vN) components located on cube sphere grid centers.
0006 %
6b4f212f67 Jean*0007 % Assume all input arrays have same shape, either original format (nc*6,nc,*)
0008 % or compact format (nc,nc*6,*), where nc is the cube face resolution.
0009 % Incoming u and v matricies are cube sphere A-grid (Grid='A') or C-grid vector
0010 % fields (defaut, Grid='C'). There may be up 2 more dimensions (likely z and
0011 % time) for a total of 4 dimensions.
0012 % Output uE,vN is returned in the shape shape (with same dimensions).
4ec37fd829 Jean*0013 % Optional maskW & maskS can be provided (for C-grid vector input case)
4582449b53 Jean*0014 % and used for cell-center interpolation
5cea1a7045 Dani*0015 %
0016 % e.g.
0017 %
0018 % >> uC=rdmds('uVeltave.0000513360');
0019 % >> vC=rdmds('vVeltave.0000513360');
0020 % >> AngleCS=rdmds('AngleCS');
0021 % >> AngleSN=rdmds('AngleSN');
0022 % >> [uE,vN] = rotate_uv2uvEN(uC,vC,AngleCS,AngleSN);
0023 %
0024 % >> uA=rdmds('uVeltaveA.0000513360');
0025 % >> vA=rdmds('vVeltaveA.0000513360');
0026 % >> AngleCS=rdmds('AngleCS');
0027 % >> AngleSN=rdmds('AngleSN');
0028 % >> [uE,vN] = rotate_uv2uvEN(uA,vA,AngleCS,AngleSN,'A');
0029
0030 % Default is a C-grid configuration.
0031 if nargin == 4, Grid = 'C'; end
4582449b53 Jean*0032 UVmsk=0; if isequal(Grid,'C') & nargin == 7, UVmsk=1; end
5cea1a7045 Dani*0033
4ec37fd829 Jean*0034 n1h=size(u,1); n2h=size(u,2);
5cea1a7045 Dani*0035 % Verify dimensions are that of cube-sphere
4ec37fd829 Jean*0036 if ~isequal(size(u),size(v)) | ...
0037 ~isequal(size(AngleCS),[n1h n2h]) | ...
0038 ~isequal(size(AngleSN),[n1h n2h])
0039 error(['Error in Input-array dimensions: ',...
5cea1a7045 Dani*0040 'u = ',mat2str(size(u)),', ',...
4ec37fd829 Jean*0041 'v = ',mat2str(size(v)),', ',...
0042 'AnCS = ',mat2str(size(AngleCS)),', ',...
0043 'AnSN = ',mat2str(size(AngleSN))]);
0044 end
0045 if n1h == n2h*6, nc=n2h; elseif n1h*6 == n2h, nc=n1h;
0046 else
0047 error(['Error in CS-dimensions: ',...
0048 'n1h,n2h = ',int2str(n1h),', ',int2str(n2h)])
0049 end
0050 dim=size(u); ncp=nc+1;
4582449b53 Jean*0051 nz=1; if length(dim) > 2, nz=prod(dim(3:end)); end
5cea1a7045 Dani*0052
4582449b53 Jean*0053 if UVmsk == 1,
0054 if length(dim) == 2,
0055 n3d=1;
a8e064706f Jean*0056 maskW=maskW(:,:,1);
0057 maskS=maskS(:,:,1);
4582449b53 Jean*0058 else
0059 n3d=dim(3);
4ec37fd829 Jean*0060 end
4582449b53 Jean*0061 dim3d=[dim(1:2) n3d];
0062 if ~isequal(size(maskW),dim3d),
0063 error(['Error in MaskW-dimensions: ',...
0064 'u = ',mat2str(size(u)),', ',...
0065 'maskW = ',mat2str(size(maskW))]);
0066 end
0067 if ~isequal(size(maskS),dim3d),
0068 error(['Error in MaskS-dimensions: ',...
0069 'v = ',mat2str(size(v)),', ',...
0070 'maskS = ',mat2str(size(maskS))]);
0071 end
0072 else
0073 dim3d=dim;
0074 end
0075
0076 %- apply mask (if provided)
0077 if UVmsk == 1,
0078 if n3d == nz,
0079 nt=1;
a8e064706f Jean*0080 u=u.*maskW;
0081 v=v.*maskS;
4582449b53 Jean*0082 else
0083 nt=dim(4);
4ec37fd829 Jean*0084 u=reshape(u,[n1h*n2h*n3d nt]);
0085 v=reshape(v,[n1h*n2h*n3d nt]);
0086 maskW=reshape(maskW,[n1h*n2h*n3d 1]);
0087 maskS=reshape(maskS,[n1h*n2h*n3d 1]);
4582449b53 Jean*0088 u=u.*(maskW*ones(1,nt));
0089 v=v.*(maskS*ones(1,nt));
0090 end
0091 end
4ec37fd829 Jean*0092 % Parse dimension information, flatten extra dimensions.
0093 u=reshape(u,[n1h n2h nz]);
0094 v=reshape(v,[n1h n2h nz]);
5cea1a7045 Dani*0095
0096 % Do simple average to put u,v at the cell center (A-grid) as needed.
0097 if isequal(Grid,'A')
4582449b53 Jean*0098 msk=ones(dim3d);
4ec37fd829 Jean*0099 u=reshape(u,[n1h*n2h nz]);
0100 v=reshape(v,[n1h*n2h nz]);
5cea1a7045 Dani*0101 elseif isequal(Grid,'C')
4582449b53 Jean*0102 [uu,vv] = split_UV_cub(u,v,1);
4ec37fd829 Jean*0103 uu=reshape(uu,[ncp nc nz 6]);
0104 vv=reshape(vv,[nc ncp nz 6]);
4582449b53 Jean*0105 if UVmsk == 1,
4ec37fd829 Jean*0106 maskW=reshape(maskW,[n1h n2h n3d]);
0107 maskS=reshape(maskS,[n1h n2h n3d]);
4582449b53 Jean*0108 [mu,mv] = split_UV_cub(maskW,maskS,0);
4ec37fd829 Jean*0109 mu=reshape(mu,[ncp nc n3d 6]);
0110 mv=reshape(mv,[nc ncp n3d 6]);
0111 um=(mu(1:nc,:,:,:)+mu(2:ncp,:,:,:))/2;
0112 vm=(mv(:,1:nc,:,:)+mv(:,2:ncp,:,:))/2;
4582449b53 Jean*0113 msk=(um+vm)/2;
0114 msk = reshape(msk,dim3d);
0115 %-----
4ec37fd829 Jean*0116 u=(uu(1:nc,:,:,:)+uu(2:ncp,:,:,:))/2;
0117 v=(vv(:,1:nc,:,:)+vv(:,2:ncp,:,:))/2;
0118 u=reshape(permute(u,[1 2 4 3]),[nc*nc*6*n3d nt]);
0119 v=reshape(permute(v,[1 2 4 3]),[nc*nc*6*n3d nt]);
0120 um=reshape(permute(um,[1 2 4 3]),[nc*nc*6*n3d 1]);
0121 vm=reshape(permute(vm,[1 2 4 3]),[nc*nc*6*n3d 1]);
0122 um(find(um==0))=1; um=1./um;
0123 vm(find(vm==0))=1; vm=1./vm;
4582449b53 Jean*0124 u=u.*(um*ones(1,nt));
0125 v=v.*(vm*ones(1,nt));
4ec37fd829 Jean*0126 if n2h == nc,
0127 u=permute(reshape(u,[nc nc 6 nz]),[1 3 2 4]);
0128 v=permute(reshape(v,[nc nc 6 nz]),[1 3 2 4]);
0129 end
4582449b53 Jean*0130 else
0131 msk=ones(dim);
4ec37fd829 Jean*0132 u=(uu(1:nc,:,:,:)+uu(2:ncp,:,:,:))/2;
0133 v=(vv(:,1:nc,:,:)+vv(:,2:ncp,:,:))/2;
0134 if n1h == nc,
0135 u=permute(u,[1 2 4 3]);
0136 v=permute(v,[1 2 4 3]);
0137 else
0138 u=permute(u,[1 4 2 3]);
0139 v=permute(v,[1 4 2 3]);
0140 end
4582449b53 Jean*0141 end
4ec37fd829 Jean*0142 u=reshape(u,[n1h*n2h nz]);
0143 v=reshape(v,[n1h*n2h nz]);
5cea1a7045 Dani*0144 else
0145 error(['Unrecognized grid type: ',Grid]);
0146 end
0147
0148 % Make rotation to find uE, vN.
4ec37fd829 Jean*0149 uE=NaN.*zeros(n1h*n2h,nz);
0150 vN=NaN.*zeros(n1h*n2h,nz);
5cea1a7045 Dani*0151 for k=1:nz;
0152 uE(:,k)=AngleCS(:).*u(:,k)-AngleSN(:).*v(:,k);
0153 vN(:,k)=AngleSN(:).*u(:,k)+AngleCS(:).*v(:,k);
0154 end
0155 uE = reshape(uE,dim);
bd73609d51 Jean*0156 vN = reshape(vN,dim);
0157 return
