| 1 |
function [fldBAR]=calc_barostream(fldU,fldV,noDiv,list_factors); |
| 2 |
%object: compute barotropic streamfunction |
| 3 |
%inputs: fldU and fldV are the fields of grid point transport |
| 4 |
%optional: noDiv (default is 1). If 1 then remove the divergent |
| 5 |
% part of the flow field first. If 0 then dont. |
| 6 |
% list_factors (default is {'dh','dz'}) |
| 7 |
%output: FLD is the streamfunction |
| 8 |
%notes: the result is converted to Sv |
| 9 |
|
| 10 |
global mygrid; |
| 11 |
|
| 12 |
if nargin<3; noDiv=1; end; |
| 13 |
if nargin<4; list_factors={'dh','dz'}; end; |
| 14 |
|
| 15 |
%0) prepare fldU/fldV (transport fields): |
| 16 |
n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1); |
| 17 |
|
| 18 |
fldU(isnan(fldU))=0; fldV(isnan(fldV))=0; |
| 19 |
|
| 20 |
dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU); |
| 21 |
if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end; |
| 22 |
facW=drf; facW(:)=1; facS=facW; |
| 23 |
for ii=1:length(list_factors); |
| 24 |
tmp1=list_factors{ii}; |
| 25 |
if strcmp(tmp1,'dh'); facW=facW.*dyg; facS=facS.*dxg; |
| 26 |
elseif strcmp(tmp1,'dz'); facW=facW.*drf; facS=facS.*drf; |
| 27 |
elseif strcmp(tmp1,'hfac'); facW=facW.*mygrid.hFacW; facS=facS.*mygrid.hFacS; |
| 28 |
elseif isempty(tmp1); 1; |
| 29 |
else; fprintf('error in calc_barostream : non supported factor\n'); return; |
| 30 |
end; |
| 31 |
end; |
| 32 |
|
| 33 |
for k4=1:n4; |
| 34 |
fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW; |
| 35 |
fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS; |
| 36 |
end; |
| 37 |
|
| 38 |
%apply mask: |
| 39 |
fldU=sum(fldU,3).*mygrid.mskW(:,:,1); |
| 40 |
fldV=sum(fldV,3).*mygrid.mskS(:,:,1); |
| 41 |
|
| 42 |
%0.1) compute streamfunction mask: |
| 43 |
if noDiv; |
| 44 |
mskU=1*~isnan(fldU); |
| 45 |
mskV=1*~isnan(fldV); |
| 46 |
[mskU,mskV]=exch_UV(mskU,mskV); |
| 47 |
mskU=abs(mskU); mskV=abs(mskV); |
| 48 |
mskBF=0*exch_T_N(mygrid.RAC); |
| 49 |
for iF=1:fldU.nFaces; |
| 50 |
tmp2=mskV{iF}; |
| 51 |
tmp3a=[ones(1,size(tmp2,2));tmp2]; |
| 52 |
tmp3b=[tmp2;ones(1,size(tmp2,2))]; |
| 53 |
mskBF{iF}=tmp3a.*tmp3b; |
| 54 |
end; |
| 55 |
end; |
| 56 |
|
| 57 |
%take out the divergent part of the flow: |
| 58 |
if noDiv; |
| 59 |
[fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV); |
| 60 |
fldU=fldU-fldUdiv; fldV=fldV-fldVdiv; |
| 61 |
end; |
| 62 |
|
| 63 |
%1) compute streamfunction face by face: |
| 64 |
[fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0; |
| 65 |
tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end; |
| 66 |
|
| 67 |
%1.1) reset one land value per face to 0: |
| 68 |
if noDiv; |
| 69 |
for iF=1:fldU.nFaces; |
| 70 |
tmpA=tmp1{iF}; |
| 71 |
tmpB=mskBF{iF}; |
| 72 |
[ii,jj]=find(tmpB==0); |
| 73 |
ii=ii(1); jj=jj(1); |
| 74 |
tmpA(:,jj)=tmpA(:,jj)-tmpB(ii,jj); |
| 75 |
for kk=1:jj-1; |
| 76 |
tmpC=tmpA(:,kk+1)-tmpA(:,kk); |
| 77 |
tmpD=fldU{iF}(:,kk); |
| 78 |
tmpE=nanmedian(tmpC+tmpD); |
| 79 |
tmpA(:,kk)=tmpA(:,kk)-tmpE; |
| 80 |
end; |
| 81 |
for kk=jj+1:size(tmpA,2); |
| 82 |
tmpC=tmpA(:,kk)-tmpA(:,kk-1); |
| 83 |
tmpD=fldU{iF}(:,kk-1); |
| 84 |
tmpE=nanmedian(tmpC+tmpD); |
| 85 |
tmpA(:,kk)=tmpA(:,kk)-tmpE; |
| 86 |
end; |
| 87 |
tmp1{iF}=tmpA; |
| 88 |
end; |
| 89 |
end; |
| 90 |
|
| 91 |
%1.2) compute divergent part of the flow on average line by line: |
| 92 |
if ~noDiv; |
| 93 |
tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1)); |
| 94 |
% to check divergence implied errors: |
| 95 |
% figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end; |
| 96 |
%subtract from streamfunction: |
| 97 |
for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end; |
| 98 |
end; |
| 99 |
|
| 100 |
bf_step1=tmp1; |
| 101 |
|
| 102 |
if fldU.nFaces==1; |
| 103 |
bf_step2=bf_step1; |
| 104 |
else; |
| 105 |
%2) match edges: |
| 106 |
%... set face number field |
| 107 |
TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end; |
| 108 |
TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points |
| 109 |
tmp1=bf_step1; |
| 110 |
for iF=1:fldU.nFaces-1; |
| 111 |
tmp2=exch_T_N(tmp1);%!!! same trick |
| 112 |
tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points |
| 113 |
TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces |
| 114 |
tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points |
| 115 |
tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points |
| 116 |
tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points |
| 117 |
tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data |
| 118 |
end; |
| 119 |
bf_step2=tmp1; |
| 120 |
end; |
| 121 |
|
| 122 |
%3) put streamfunction at cell center |
| 123 |
tmp1=bf_step2; |
| 124 |
tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2; |
| 125 |
tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end; |
| 126 |
bf_step3=tmp2; |
| 127 |
%4) set 0 on fist land point: |
| 128 |
tmp1=convert2vector(bf_step3); |
| 129 |
tmp2=convert2vector(mygrid.mskC(:,:,1)); |
| 130 |
tmp2=find(isnan(tmp2)&~isnan(tmp1)); |
| 131 |
% |
| 132 |
%tmp2=median(tmp1(tmp2));%the original method |
| 133 |
%tmp2=tmp1(tmp2(1));%a point in Antarctica (in LLC90 at least) |
| 134 |
%closest land point closest to Boston: |
| 135 |
tmp_lon=convert2vector(mygrid.XC); tmp_lon=tmp_lon(tmp2); |
| 136 |
tmp_lat=convert2vector(mygrid.YC); tmp_lat=tmp_lat(tmp2); |
| 137 |
tmp_dis=(tmp_lat-42.3601).^2+(tmp_lon-71.0589).^2; |
| 138 |
tmp2=tmp1(tmp2(find(tmp_dis==min(tmp_dis)))); |
| 139 |
% |
| 140 |
bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1); |
| 141 |
|
| 142 |
%5) return the result: |
| 143 |
fldBAR=bf_step4; |
| 144 |
|
| 145 |
%convert to Sv and change sign: |
| 146 |
fldBAR=1e-6*fldBAR; |
| 147 |
|
| 148 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph