Warning, /utils/matlab/cs_grid/calcBolusPsiCube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 0ff2855e on 2006-08-18 00:51:40 UTC
0ff2855e12 Davi*0001 function [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile);
0002
0003 % [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile);
0004 %
0005 % Compute bolus streamfunction from GM scheme
0006 %
0007 % Input arguments:
0008 % The incoming field data (d) and grid data (g) must be in a structured
0009 % array format (which is the format that comes from rdmnc):
0010 % d [Field data] Kwx,Kwy
0011 % g [Grid data ] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
0012 % GMform [string] GM form 'Skew' or 'Advc'
0013 % blkFile [file name] Broken line file
0014 % Output arguments:
0015 % PsiB : bolus streamfunction at interface level (in Sv)
0016 % ylat : meridional coordinate of PsiB
0017 %
0018 % Comments:
0019 % -For Skew-flux form:
0020 % PsiB computed from Kwx & Kwy divided by 2.
0021 % first average Kwx and Kwy at u- and v-points:
0022 % psiX=(rAc*Kwx)_i / (dXc*dYg) ; psiY=(rAc*Kwy)_j / dYc ;
0023 % and then "zonally" average along broken lines
0024 % -For Advective form:
0025 % PsiB computed from PsiX and PsiY
0026 % just need to "zonally" average along broken lines
0027 %
0028 %---------------------------------------------------------------------
0029
0030
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 % Prepare grid stuff %
0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034
0035 nc = size(g.XC,2);
0036 nr = length(g.drF);
0037 nt = size(d.GM_Kwx,4);
0038
0039 %--- areas :
0040 ra = g.rA;
0041 dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);
0042 dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);
0043 dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);
0044 dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]);
0045
0046 rAu=dxc.*dyg;
0047 rAv=dyc.*dxg;
0048
0049 %--- masks :
0050 hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0051 hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0052 mskw=ceil(hw); mskw=min(1,mskw);
0053 msks=ceil(hs); msks=min(1,msks);
0054
0055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0056 % Read/Prepare GM fields %
0057 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0058 psiX_all = zeros(6*nc*nc,nr,nt);
0059 psiY_all = zeros(6*nc*nc,nr,nt);
0060
0061 switch GMform
0062 case 'Skew'
0063
0064 kwx_all = 0.5*d.GM_Kwx;
0065 kwy_all = 0.5*d.GM_Kwy;
0066
0067 for it = 1:nt
0068 kwx = kwx_all(:,:,:,it);
0069 kwy = kwy_all(:,:,:,it);
0070
0071 %-- K*ra + add 1 overlap :
0072 kwx = repmat(ra,[1 1 nr]).*kwx;
0073 kwy = repmat(ra,[1 1 nr]).*kwy;
0074 v6X = split_C_cub(kwx,1);
0075 v6Y = split_C_cub(kwy,1);
0076 k6x = v6X(:,[2:nc+1],:,:);
0077 k6y = v6Y([2:nc+1],:,:,:);
0078
0079 %-----------------
0080 clear v6X v6Y
0081 v6X = (k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:))/2;
0082 v6Y = (k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:))/2;
0083
0084 psiX = zeros(6*nc,nc,nr);
0085 psiY = zeros(6*nc,nc,nr);
0086
0087 for n = 1:6
0088 is = 1+nc*(n-1);ie=nc*n;
0089 psiX([is:ie],[1:nc],[1:nr]) = v6X(:,:,:,n);
0090 psiY([is:ie],[1:nc],[1:nr]) = v6Y(:,:,:,n);
0091 end
0092
0093 psiX = reshape(psiX,6*nc*nc,nr);
0094 psiY = reshape(psiY,6*nc*nc,nr);
0095
0096 psiX_all(:,:,it) = mskw .* psiX ./ repmat(rAu,[1,nr]);
0097 psiY_all(:,:,it) = msks .* psiY ./ repmat(rAv,[1,nr]);
0098
0099 end
0100
0101 case 'Advc'
0102
0103 psiX_all = reshape(d.GM_PsiX(1:6*nc,:,:,1:nt),6*nc*nc,nr,nt);
0104 psiY_all = reshape(d.GM_PsiY(:,1:nc,:,1:nt) ,6*nc*nc,nr,nt);
0105
0106 otherwise
0107 disp(['C est Portnawak: GMform should be Skew or Advc: ',GMform])
0108 end
0109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0110 % Zonally integrate along broken lines %
0111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0112
0113 load(blkFile);
0114 ydim = length(bkl_Ylat);
0115 ylat = [-90,bkl_Ylat,90];
0116 ufac= rem(bkl_Flg,2);
0117 vfac= fix(bkl_Flg/2);
0118
0119 PsiB= zeros(ydim+2,nr+1,nt);
0120
0121 for it = 1:nt
0122 for k = 1:nr
0123 psixt=dyg.*psiX_all(:,k,it);
0124 psiyt=dxg.*psiY_all(:,k,it);
0125 for jl = 1:ydim
0126 ie = bkl_Npts(jl);
0127 PsiB(jl+1,k,it) = sum( ufac(1:ie,jl).*psixt(bkl_IJuv(1:ie,jl)) ...
0128 + vfac(1:ie,jl).*psiyt(bkl_IJuv(1:ie,jl)) );
0129 end
0130 end
0131 end
0132 PsiB = 1e-6*PsiB;
