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

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

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

Revision 1.7 - (hide annotations) (download)
Mon Sep 24 23:05:20 2012 UTC (12 years, 10 months ago) by gforget
Branch: MAIN
Changes since 1.6: +28 -9 lines 
- calc_barostream.m : add list_factors argument to accomodate e.g. tracer class transports.
- calc_overturn.m : add doFlip,list_factors args to accomodate e.g. tracer class transports.
- calc_barostream.m : introduce list_factors to allow use for e.g. tracer class transports.
- layers_remap.m (new) : remap variables (e.g. transports) from depth to tracer classes.
                       Uses regrid_dblres.m and regrid_sum.m and mimics pkg/layers.
- regrid_dblres.m (new) : double the resolution (only along 3rd dimension for now)
                       for a variable P (extensive or intensive) a number of times.
- regrid_sum.m (new) : add 3rd dimension elements of extensive variable P,
                       according to values of a tracer field collocated with P,
                       to the tracer grid defined by trGrid (1D vector)
1 gforget 1.7 function [fldBAR]=calc_barostream(fldU,fldV,noDiv,list_factors);
2 gforget 1.5 %object: compute barotropic streamfunction
3     %inputs: fldU and fldV are the fields of grid point transport
4 gforget 1.7 %optional: noDiv (default is 1). If 1 then remove the divergent
5     % part of the flow field first. If 0 then dont.
6     % list_factors (default is {'dh','dz'})
7 gforget 1.5 %output: FLD is the streamfunction
8 gforget 1.7 %notes: the result is converted to Sv
9 gforget 1.1
10     global mygrid;
11    
12 gforget 1.7 if nargin<3; noDiv=1; end;
13     if nargin<4; list_factors={'dh','dz'}; end;
14    
15 gforget 1.2 %0) prepare fldU/fldV (transport fields):
16 gforget 1.1 n3=max(size(fldU.f1,3),1); n4=max(size(fldV.f1,4),1);
17    
18     fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
19    
20 gforget 1.7
21     dxg=mk3D(mygrid.DXG,fldU); dyg=mk3D(mygrid.DYG,fldU);
22     if size(fldU.f1,3)==length(mygrid.DRF); drf=mk3D(mygrid.DRF,fldU); else; drf=fldU; drf(:)=1; end;
23     facW=drf; facW(:)=1; facS=facW;
24     for ii=1:length(list_factors);
25     tmp1=list_factors{ii};
26     if strcmp(tmp1,'dh'); facW=facW.*dyg; facS=facS.*dxg;
27     elseif strcmp(tmp1,'dz'); facW=facW.*drf; facS=facS.*drf;
28     elseif strcmp(tmp1,'hfac'); facW=facW.*mygrid.hFacW; facS=facS.*mygrid.hFacS;
29     elseif isempty(tmp1); 1;
30     else; fprintf('error in calc_UV_div: non supported factor\n'); return;
31     end;
32     end;
33    
34 gforget 1.1 for k4=1:n4;
35 gforget 1.7 fldU(:,:,:,k4)=fldU(:,:,:,k4).*facW;
36     fldV(:,:,:,k4)=fldV(:,:,:,k4).*facS;
37 gforget 1.1 end;
38    
39 gforget 1.3 %apply mask:
40     fldU=sum(fldU,3).*mygrid.mskW(:,:,1);
41     fldV=sum(fldV,3).*mygrid.mskS(:,:,1);
42     %take out thedivergent part of the flow:
43 gforget 1.7 if noDiv;
44 gforget 1.3 [fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV);
45     fldU=fldU-fldUdiv; fldV=fldV-fldVdiv;
46     end;
47 gforget 1.1
48 gforget 1.2 %1) compute streamfunction face by face:
49     [fldU,fldV]=exch_UV(fldU,fldV); fldU(isnan(fldU))=0; fldV(isnan(fldV))=0;
50     tmp1=cumsum(fldV,1); for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=[zeros(1,size(tmp2,2));tmp2]; end;
51     tmp2=diff(tmp1,1,2)+fldU; tmp3=cumsum(mean(tmp2,1));
52 gforget 1.3 %to check divergence implied errors:
53     %figure; for iF=1:fldU.nFaces; subplot(3,2,iF); plot(std(tmp2{iF},0,1)); end;
54 gforget 1.2 for iF=1:fldU.nFaces; tmp2=tmp1{iF}; tmp1{iF}=tmp2-ones(size(tmp2,1),1)*[0 tmp3{iF}]; end;
55     bf_step1=tmp1;
56    
57 gforget 1.3 if fldU.nFaces==1;
58 gforget 1.4 bf_step2=bf_step1;
59 gforget 1.3 else;
60 gforget 1.2 %2) match edges:
61     %... set face number field
62     TMP1=tmp1; for iF=1:TMP1.nFaces; TMP1{iF}(:)=iF; end;
63     TMP2=exch_T_N(TMP1);%!!! this is a trick, since TMP1 is (n+1 X n+1) and loc. at vorticity points
64     tmp1=bf_step1;
65     for iF=1:fldU.nFaces-1;
66     tmp2=exch_T_N(tmp1);%!!! same trick
67     tmp3=tmp2{iF+1}; tmp3(3:end-2,3:end-2)=NaN;%mask out interior points
68     TMP3=TMP2{iF+1}; tmp3(find(TMP3>iF+1))=NaN;%mask out edges points coming from unadjusted faces
69     tmp3(:,1)=tmp3(:,1)-tmp3(:,2); tmp3(:,end)=tmp3(:,end)-tmp3(:,end-1);%compare edge points
70     tmp3(1,:)=tmp3(1,:)-tmp3(2,:); tmp3(end,:)=tmp3(end,:)-tmp3(end-1,:);%compare edge points
71     tmp3(2:end-1,2:end-1)=NaN;%mask out remaining interior points
72     tmp1{iF+1}=tmp1{iF+1}+nanmedian(tmp3(:));%adjust the face data
73 gforget 1.1 end;
74 gforget 1.2 bf_step2=tmp1;
75 gforget 1.4 end;
76 gforget 1.2 %3) put streamfunction at cell center
77     tmp1=bf_step2;
78     tmp2=tmp1; for iF=1:tmp1.nFaces; tmp3=tmp2{iF}; tmp3=(tmp3(:,1:end-1)+tmp3(:,2:end))/2;
79     tmp3=(tmp3(1:end-1,:)+tmp3(2:end,:))/2; tmp2{iF}=tmp3; end;
80     bf_step3=tmp2;
81     %4) set 0 on land on average:
82     tmp1=convert2array(bf_step3);
83     tmp2=convert2array(mygrid.mskC(:,:,1)); tmp2=find(isnan(tmp2)&~isnan(tmp1));
84     tmp2=median(tmp1(tmp2)); if isnan(tmp2); tmp2=nanmedian(tmp1(:)); end;
85     bf_step4=(bf_step3-tmp2).*mygrid.mskC(:,:,1);
86 gforget 1.1
87 gforget 1.2 %5) return the result:
88     fldBAR=bf_step4;
89 gforget 1.1
90 gforget 1.7 %convert to Sv and change sign:
91     fldBAR=1e-6*fldBAR;
92    
  ViewVC Help
Powered by ViewVC 1.1.22