Warning, /utils/matlab/cs_grid/calc_vort_cs.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
896b20e91e Jean*0001 function [vort,z6t]=calc_vort_cs(u3d,v3d,dxC,dyC,rAz);
0002 % [vort,z6t]=calc_vort_cs(u3d,v3d,dxC,dyC,rAz);
0003 % compute vorticity (3rd component) on CS-grid.
6b4f212f67 Jean*0004 % assume all input (except rAz) have same shape, either original format (nc*6,nc,*)
0005 % or compact format (nc,nc*6,*); and rAz can either have same shape as others
0006 % or can include 2 missing corner in a long vector shape (nc*nc*6+2)
896b20e91e Jean*0007 % output is provided with 2 shapes:
6b4f212f67 Jean*0008 % 1) vort(nc*nc*6+2,*) = long-vector (compact) form;
0009 % 2) z6t(nc+1,nc+1,*,6) = face splitted.
896b20e91e Jean*0010 %
6b4f212f67 Jean*0011 % Written by jmc@mit.edu, 2005.
0012
0013 dims=size(u3d);
0014 n1h=dims(1); n2h=dims(2);
0015 %-- Verify input array dimensions
0016 if n1h == n2h*6, nc=n2h; elseif n1h*6 == n2h, nc=n1h;
0017 else
0018 error(['Error in CS-dimensions: ',...
0019 'n1h,n2h = ',int2str(n1h),', ',int2str(n2h)])
0020 end
0021 ncp=nc+1; n2p=nc+2; nPg=nc*nc*6 ;
0022 if ~isequal(size(u3d),size(v3d)) | ...
0023 ~isequal(size(dxC),[n1h n2h]) | ...
0024 ~isequal(size(dyC),[n1h n2h])
0025 error(['Error in Input-array dimensions: ',...
0026 'u3d = ',mat2str(size(u3d)),', ',...
0027 'v3d = ',mat2str(size(v3d)),', ',...
0028 'dxC = ',mat2str(size(dxC)),', ',...
0029 'dyC = ',mat2str(size(dyC))]);
0030 end
0031 if ~isequal(size(rAz),[n1h n2h]) & ~isequal(size(rAz),[nPg+2 1]),
0032 error(['Error in rAz array dimensions: ',mat2str(size(rAz))]);
0033 end
0034
896b20e91e Jean*0035 if length(dims) == 2, nr=1; else nr=prod(dims(3:end)); end
0036
0037 siZ=prod(size(rAz));
6b4f212f67 Jean*0038 if size(rAz,1) == nPg+2,
896b20e91e Jean*0039 rAz=reshape(rAz,[nPg+2 1]);
0040 aZ=split_Z_cub(rAz);
6b4f212f67 Jean*0041 else
0042 if n2h == nc,
0043 rAz=permute(reshape(rAz,[nc 6 nc]),[1 3 2]);
0044 end
0045 rAz=reshape(rAz,[nPg 1]);
896b20e91e Jean*0046 aZc=zeros(nPg+2,1); aZc(1:nPg,:)=rAz; aZc(nPg+1)=rAz(1); aZc(nPg+2)=rAz(1);
0047 aZ=split_Z_cub(aZc);
0048 end
0049
6b4f212f67 Jean*0050 u3d=reshape(u3d,[n1h n2h nr]);
0051 v3d=reshape(v3d,[n1h n2h nr]);
896b20e91e Jean*0052
0053 [u6t,v6t]=split_UV_cub(u3d,v3d,1,2);
0054 [d6x,d6y]=split_UV_cub(dxC,dyC,0,2);
6b4f212f67 Jean*0055 if nr > 1,
896b20e91e Jean*0056 u6t=permute(u6t,[1 2 4 3]);
0057 v6t=permute(v6t,[1 2 4 3]);
6b4f212f67 Jean*0058 end
896b20e91e Jean*0059 z6t=zeros(ncp,ncp,6,nr);
0060
0061 for k=1:nr,
0062 vv1=d6x.*u6t(:,:,:,k);
0063 vv2=d6y.*v6t(:,:,:,k);
0064 z6t(:,:,:,k)=vv2([2:n2p],:,:)-vv2([1:ncp],:,:);
0065 z6t(:,:,:,k)=z6t(:,:,:,k)-(vv1(:,[2:n2p],:)-vv1(:,[1:ncp],:));
0066 %- corner: the quick way:
0067 % z6t(1,1,:,k) = vv2(2,1,:) -vv2(1,1,:) -vv1(1,2,:);
0068 % z6t(1,ncp,:,k) = vv2(2,ncp,:)-vv2(1,ncp,:)+vv1(1,ncp,:);
0069 % z6t(ncp,1,:,k) = vv2(n2p,1,:)-vv2(ncp,1,:)-vv1(ncp,2,:);
0070 % z6t(ncp,ncp,:,k)=vv2(n2p,ncp,:)-vv2(ncp,ncp,:)+vv1(ncp,ncp,:);
6b4f212f67 Jean*0071 %- corner: add the 3 terms always in the same order
896b20e91e Jean*0072 % to get the same truncation on the 3 faces
0073 for n=1:3,
0074 f=2*n-1; %- odd face number
0075 vc=-vv2(1,1,f); %- S-W corner
e61bf3af22 Jean*0076 z6t(1,1,f,k) = (vv2(2,1,f)-vv1(1,2,f))+vc;
896b20e91e Jean*0077 vc=+vv2(n2p,1,f); %- S-E corner
e61bf3af22 Jean*0078 z6t(ncp,1,f,k) = (vc-vv1(ncp,2,f))-vv2(ncp,1,f);
896b20e91e Jean*0079 vc=+vv2(n2p,ncp,f); %- N-E corner
e61bf3af22 Jean*0080 z6t(ncp,ncp,f,k)=(vc-vv2(ncp,ncp,f))+vv1(ncp,ncp,f);
896b20e91e Jean*0081 vc=-vv2(1,ncp,f); %- N-W corner
e61bf3af22 Jean*0082 vc3=[vc vv2(2,ncp,f) vv1(1,ncp,f) vc vv2(2,ncp,f)];
896b20e91e Jean*0083 z6t(1,ncp,f,k) = (vc3(n+2)+vc3(n+1))+vc3(n);
0084 f=2*n; %- even face number
0085 vc=-vv2(1,1,f); %- S-W corner
e61bf3af22 Jean*0086 z6t(1,1,f,k) = (vv2(2,1,f)-vv1(1,2,f))+vc;
896b20e91e Jean*0087 vc=+vv2(n2p,1,f); %- S-E corner
0088 vc3=[-vv1(ncp,2,f) -vv2(ncp,1,f) vc -vv1(ncp,2,f) -vv2(ncp,1,f)];
0089 z6t(ncp,1,f,k) = (vc3(n)+vc3(n+1))+vc3(n+2);
0090 vc=+vv2(n2p,ncp,f); %- N-E corner
e61bf3af22 Jean*0091 z6t(ncp,ncp,f,k)=(vv1(ncp,ncp,f)+vc)-vv2(ncp,ncp,f);
896b20e91e Jean*0092 vc=-vv2(1,ncp,f); %- N-W corner
e61bf3af22 Jean*0093 z6t(1,ncp,f,k) = (vv2(2,ncp,f)+vc)+vv1(1,ncp,f);
896b20e91e Jean*0094 end
0095 %- divide by rAz:
0096 z6t(:,:,:,k)=z6t(:,:,:,k)./aZ;
0097 end
0098
6b4f212f67 Jean*0099 %- put in compressed form:
896b20e91e Jean*0100 vort=zeros(nPg+2,nr);
0101 % extract the interior
6b4f212f67 Jean*0102 vort([1:nPg],:)=reshape(z6t(1:nc,1:nc,:,:),[nPg nr]);
896b20e91e Jean*0103 % put back the 2 missing corners (N.W corner of 1rst face & S.E corner of 2nd face)
0104 vort(nPg+1,:)=z6t(1,ncp,1,:);
0105 vort(nPg+2,:)=z6t(ncp,1,2,:);
6b4f212f67 Jean*0106
896b20e91e Jean*0107 %- back into final shape:
0108 z6t=permute(z6t,[1 2 4 3]);
6b4f212f67 Jean*0109 if length(dims) == 2,
896b20e91e Jean*0110 vort=squeeze(vort);
0111 z6t=squeeze(z6t);
0112 else
0113 vort=reshape(vort,[nPg+2 dims(3:end)]);
0114 z6t=reshape(z6t,[ncp ncp dims(3:end) 6]);
0115 end
0116
6b4f212f67 Jean*0117 %-----------------
896b20e91e Jean*0118 return
