/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_bolus.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_bolus.m

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


Revision 1.2 - (hide annotations) (download)
Fri Apr 17 13:06:16 2015 UTC (10 years, 3 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65r, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
Changes since 1.1: +4 -0 lines
- fix problematic cases

1 gforget 1.1 function [bolusU,bolusV,bolusW]=calc_bolus(GM_PsiX,GM_PsiY);
2     %object: compute bolus velocty field (bolusU,bolusV,bolusW)
3     % from gm streamfunction (GM_PsiX,GM_PsiY).
4     %input: GM_PsiX,GM_PsiY
5     %output: bolusU,bolusV,bolusW
6    
7     gcmfaces_global;
8     nr=length(mygrid.RC);
9    
10 gforget 1.2 %replace NaNs with 0s:
11     GM_PsiX(isnan(GM_PsiX))=0;
12     GM_PsiY(isnan(GM_PsiY))=0;
13    
14 gforget 1.1 %compute bolus velocity:
15     GM_PsiX(:,:,nr+1)=0*GM_PsiX(:,:,end);
16     GM_PsiY(:,:,nr+1)=0*GM_PsiY(:,:,end);
17     bolusU=0*mygrid.hFacW;
18     bolusV=0*mygrid.hFacS;
19     for k=1:nr;
20     bolusU(:,:,k)=(GM_PsiX(:,:,k+1)-GM_PsiX(:,:,k))/mygrid.DRF(k);
21     bolusV(:,:,k)=(GM_PsiY(:,:,k+1)-GM_PsiY(:,:,k))/mygrid.DRF(k);
22     end;
23     bolusU=bolusU.*(~isnan(mygrid.mskW));
24     bolusV=bolusV.*(~isnan(mygrid.mskS));
25    
26     %and its vertical part
27     % (seems correct, leading to 0 divergence)
28     tmp_x=GM_PsiX(:,:,1:nr).*repmat(mygrid.DYG,[1 1 nr]);
29     tmp_y=GM_PsiY(:,:,1:nr).*repmat(mygrid.DXG,[1 1 nr]);
30     [tmp_x,tmp_y]=exch_UV(tmp_x,tmp_y);
31     tmp_a=repmat(mygrid.RAC,[1 1 nr]);
32     tmp_w=0*tmp_x;
33     for iFace=1:mygrid.nFaces;
34     tmp_w{iFace}=( tmp_x{iFace}(2:end,:,:)-tmp_x{iFace}(1:end-1,:,:) )+...
35     ( tmp_y{iFace}(:,2:end,:)-tmp_y{iFace}(:,1:end-1,:) );
36     tmp_w{iFace}=tmp_w{iFace}./tmp_a{iFace};
37     end;
38     bolusW=tmp_w.*(~isnan(mygrid.mskC));

  ViewVC Help
Powered by ViewVC 1.1.22