--- MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m 2013/03/05 18:08:05 1.1 +++ MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m 2018/04/19 00:03:59 1.3 @@ -1,4 +1,4 @@ -function [ub,vb]=calcBolusVelCube(d,g,GMform); +function [ub,vb,wb]=calcBolusVelCube(d,g,GMform); % [ub,vb] = calcBolusVelCube(d,g,GMform); % @@ -44,22 +44,23 @@ ra = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]); rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]); rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]); -%dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]); -%dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]); -%dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); -%dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); - -%rAw=dxc.*dyg; -%rAs=dyc.*dxg; +dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]); +dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]); +dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); +dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); %--- recip_hFac & mask : mw=ceil(hw); mw=min(1,mw); ms=ceil(hs); ms=min(1,ms); -%hw(find(hw==0))=Inf; -%hs(find(hs==0))=Inf; -%hw_recip=1./hw; %hw_recip(find(hw==0))=0; -%hs_recip=1./hs; %hs_recip(find(hs==0))=0; +hw(find(hw==0))=Inf; +hs(find(hs==0))=Inf; +hw_recip=1./hw; %hw_recip(find(hw==0))=0; +hs_recip=1./hs; %hs_recip(find(hs==0))=0; + +ub_all = zeros(6*nc,nc,nr,nt); +vb_all = zeros(6*nc,nc,nr,nt); +wb_all = zeros(6*nc,nc,nr,nt); switch GMform @@ -117,9 +118,42 @@ ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); % vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); - + ub_all(:,:,:,it) = ub; vb_all(:,:,:,it) = vb; + +%%%%%%%%%%%%% + [u6t,v6t] =split_UV_cub(ub,vb,0,1); + [hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1); + [dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1); + + %F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); + %F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); + F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); + F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); + + Hdiv = zeros(nc,nc,nr,6); + Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ... + + F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:); + for k=1:nr + Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k); + end + + %psiX = zeros(6*nc,nc,nr); + %for n = 1:6 + % is = 1+nc*(n-1);ie=nc*n; + % psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n); + %end + psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr); + %wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]); + + wb = zeros(6*nc,nc,nr+1); + for k=nr:-1:1 + wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1); + end + + wb_all(:,:,:,it) = wb(:,:,1:30); +%%%%%%%%%%%% end %%%%%%% Advective form case @@ -160,3 +194,4 @@ ub = ub_all; vb = vb_all; +wb = wb_all;