| 1 |
gforget |
1.3 |
function [fldBAR]=calc_barostream(fldU,fldV,varargin); |
| 2 |
gforget |
1.1 |
|
| 3 |
|
|
%by assumption: grid_load is done |
| 4 |
|
|
global mygrid; |
| 5 |
|
|
|
| 6 |
gforget |
1.2 |
%0) prepare fldU/fldV (transport fields): |
| 7 |
gforget |
1.1 |
n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1); |
| 8 |
|
|
|
| 9 |
|
|
fldU(isnan(fldU))=0; fldV(isnan(fldV))=0; |
| 10 |
|
|
|
| 11 |
|
|
dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU); drf=mk3D(mygrid.DRF,fldU); |
| 12 |
|
|
for k4=1:n4; |
| 13 |
|
|
fldU(:,:,:,k4)=fldU(:,:,:,k4).*dyg.*drf.*mygrid.hFacW*1e-6; |
| 14 |
|
|
fldV(:,:,:,k4)=fldV(:,:,:,k4).*dxg.*drf.*mygrid.hFacS*1e-6; |
| 15 |
|
|
end; |
| 16 |
|
|
|
| 17 |
gforget |
1.3 |
%apply mask: |
| 18 |
|
|
fldU=sum(fldU,3).*mygrid.mskW(:,:,1); |
| 19 |
|
|
fldV=sum(fldV,3).*mygrid.mskS(:,:,1); |
| 20 |
|
|
%take out thedivergent part of the flow: |
| 21 |
|
|
if nargin==3; doRemoveDivergentPart=varargin{1}; else; doRemoveDivergentPart=1; end; |
| 22 |
|
|
if doRemoveDivergentPart; |
| 23 |
|
|
[fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV); |
| 24 |
|
|
fldU=fldU-fldUdiv; fldV=fldV-fldVdiv; |
| 25 |
|
|
end; |
| 26 |
gforget |
1.1 |
|
| 27 |
gforget |
1.2 |
%1) compute streamfunction face by face: |
| 28 |
|
|
[fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0; |
| 29 |
|
|
tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end; |
| 30 |
|
|
tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1)); |
| 31 |
gforget |
1.3 |
%to check divergence implied errors: |
| 32 |
|
|
%figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end; |
| 33 |
gforget |
1.2 |
for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end; |
| 34 |
|
|
bf_step1=tmp1; |
| 35 |
|
|
|
| 36 |
gforget |
1.3 |
if fldU.nFaces==1; |
| 37 |
gforget |
1.4 |
bf_step2=bf_step1; |
| 38 |
gforget |
1.3 |
else; |
| 39 |
gforget |
1.2 |
%2) match edges: |
| 40 |
|
|
%... set face number field |
| 41 |
|
|
TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end; |
| 42 |
|
|
TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points |
| 43 |
|
|
tmp1=bf_step1; |
| 44 |
|
|
for iF=1:fldU.nFaces-1; |
| 45 |
|
|
tmp2=exch_T_N(tmp1);%!!! same trick |
| 46 |
|
|
tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points |
| 47 |
|
|
TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces |
| 48 |
|
|
tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points |
| 49 |
|
|
tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points |
| 50 |
|
|
tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points |
| 51 |
|
|
tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data |
| 52 |
gforget |
1.1 |
end; |
| 53 |
gforget |
1.2 |
bf_step2=tmp1; |
| 54 |
gforget |
1.4 |
end; |
| 55 |
gforget |
1.2 |
%3) put streamfunction at cell center |
| 56 |
|
|
tmp1=bf_step2; |
| 57 |
|
|
tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2; |
| 58 |
|
|
tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end; |
| 59 |
|
|
bf_step3=tmp2; |
| 60 |
|
|
%4) set 0 on land on average: |
| 61 |
|
|
tmp1=convert2array(bf_step3); |
| 62 |
|
|
tmp2=convert2array(mygrid.mskC(:,:,1)); tmp2=find(isnan(tmp2)&~isnan(tmp1)); |
| 63 |
|
|
tmp2=median(tmp1(tmp2)); if isnan(tmp2); tmp2=nanmedian(tmp1(:)); end; |
| 64 |
|
|
bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1); |
| 65 |
gforget |
1.1 |
|
| 66 |
gforget |
1.2 |
%5) return the result: |
| 67 |
|
|
fldBAR=bf_step4; |
| 68 |
gforget |
1.1 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph