/[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.2 by dfer, Wed Mar 7 22:01:18 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]);
51    
52  %rAw=dxc.*dyg;  %rAw=dxc.*dyg;
53  %rAs=dyc.*dxg;  %rAs=dyc.*dxg;
# Line 56  rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc* Line 56  rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*
56  mw=ceil(hw); mw=min(1,mw);  mw=ceil(hw); mw=min(1,mw);
57  ms=ceil(hs); ms=min(1,ms);  ms=ceil(hs); ms=min(1,ms);
58    
59  %hw(find(hw==0))=Inf;  hw(find(hw==0))=Inf;
60  %hs(find(hs==0))=Inf;  hs(find(hs==0))=Inf;
61  %hw_recip=1./hw; %hw_recip(find(hw==0))=0;  hw_recip=1./hw; %hw_recip(find(hw==0))=0;
62  %hs_recip=1./hs; %hs_recip(find(hs==0))=0;  hs_recip=1./hs; %hs_recip(find(hs==0))=0;
63    
64    ub_all = zeros(6*nc,nc,nr,nt);
65    vb_all = zeros(6*nc,nc,nr,nt);
66    wb_all = zeros(6*nc,nc,nr,nt);
67    
68  switch GMform  switch GMform
69    
# Line 117  for it = 1:nt Line 121  for it = 1:nt
121      ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);      ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]);
122  %    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);  %    vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]);
123      vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);      vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]);
124        
125      ub_all(:,:,:,it) = ub;      ub_all(:,:,:,it) = ub;
126      vb_all(:,:,:,it) = vb;      vb_all(:,:,:,it) = vb;
127    
128    %%%%%%%%%%%%%
129        [u6t,v6t]  =split_UV_cub(ub,vb,0,1);
130        [hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1);
131        [dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1);
132    
133        %F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
134        %F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);
135        F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]);
136        F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]);
137    
138        Hdiv = zeros(nc,nc,nr,6);
139        Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ...
140             + F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:);
141        for k=1:nr
142          Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k);
143        end
144    
145        %psiX = zeros(6*nc,nc,nr);
146        %for n = 1:6
147        %  is = 1+nc*(n-1);ie=nc*n;
148        %  psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n);
149        %end
150        psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr);
151        %wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]);
152    
153        wb = zeros(6*nc,nc,nr+1);
154        for k=nr:-1:1
155         wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1);
156        end
157    
158        wb_all(:,:,:,it) = wb(:,:,1:30);
159    %%%%%%%%%%%%
160  end  end
161    
162  %%%%%%% Advective form case  %%%%%%% Advective form case
# Line 160  end Line 197  end
197    
198  ub = ub_all;  ub = ub_all;
199  vb = vb_all;  vb = vb_all;
200    wb = wb_all;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22