Diff of /MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m

Parent Directory | Revision Log | Revision Graph | Patch

--- MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m	2013/03/05 18:08:05	1.1
+++ MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m	2018/03/07 22:01:18	1.2
@@ -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,10 +44,10 @@
 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]);
+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;
@@ -56,10 +56,14 @@
 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 +121,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 +197,4 @@
 
 ub = ub_all;
 vb = vb_all;
+wb = wb_all;