/[MITgcm]/MITgcm_contrib/dfer/matlab_stuff/calcBolusVelCube.m
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by dfer, Tue Mar 5 18:08:05 2013 UTC revision 1.3 by dfer, Thu Apr 19 00:03:59 2018 UTC
# Line 1  Line 1 
1  function [ub,vb]=calcBolusVelCube(d,g,GMform);  function [ub,vb,wb]=calcBolusVelCube(d,g,GMform);
2    
3  % [ub,vb] = calcBolusVelCube(d,g,GMform);  % [ub,vb] = calcBolusVelCube(d,g,GMform);
4  %  %
# Line 44  hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[ Line 44  hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[
44  ra  = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]);  ra  = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]);
45  rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]);  rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]);
46  rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]);  rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]);
47  %dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);  dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);
48  %dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);  dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);
49  %dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);  dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);
50  %dyg = reshape(g.dyG(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;  
51    
52  %--- recip_hFac & mask :  %--- recip_hFac & mask :
53  mw=ceil(hw); mw=min(1,mw);  mw=ceil(hw); mw=min(1,mw);
54  ms=ceil(hs); ms=min(1,ms);  ms=ceil(hs); ms=min(1,ms);
55    
56  %hw(find(hw==0))=Inf;  hw(find(hw==0))=Inf;
57  %hs(find(hs==0))=Inf;  hs(find(hs==0))=Inf;
58  %hw_recip=1./hw; %hw_recip(find(hw==0))=0;  hw_recip=1./hw; %hw_recip(find(hw==0))=0;
59  %hs_recip=1./hs; %hs_recip(find(hs==0))=0;  hs_recip=1./hs; %hs_recip(find(hs==0))=0;
60    
61    ub_all = zeros(6*nc,nc,nr,nt);
62    vb_all = zeros(6*nc,nc,nr,nt);
63    wb_all = zeros(6*nc,nc,nr,nt);
64    
65  switch GMform  switch GMform
66    
# Line 117  for it = 1:nt Line 118  for it = 1:nt
118      ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);      ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);
119  %    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);  %    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);
120      vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);      vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);
121        
122      ub_all(:,:,:,it) = ub;      ub_all(:,:,:,it) = ub;
123      vb_all(:,:,:,it) = vb;      vb_all(:,:,:,it) = vb;
124    
125    %%%%%%%%%%%%%
126        [u6t,v6t]  =split_UV_cub(ub,vb,0,1);
127        [hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1);
128        [dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1);
129    
130        %F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
131        %F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);
132        F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
133        F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);
134    
135        Hdiv = zeros(nc,nc,nr,6);
136        Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ...
137             + F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:);
138        for k=1:nr
139          Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k);
140        end
141    
142        %psiX = zeros(6*nc,nc,nr);
143        %for n = 1:6
144        %  is = 1+nc*(n-1);ie=nc*n;
145        %  psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n);
146        %end
147        psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr);
148        %wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]);
149    
150        wb = zeros(6*nc,nc,nr+1);
151        for k=nr:-1:1
152         wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1);
153        end
154    
155        wb_all(:,:,:,it) = wb(:,:,1:30);
156    %%%%%%%%%%%%
157  end  end
158    
159  %%%%%%% Advective form case  %%%%%%% Advective form case
# Line 160  end Line 194  end
194    
195  ub = ub_all;  ub = ub_all;
196  vb = vb_all;  vb = vb_all;
197    wb = wb_all;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22