Warning, /utils/matlab/cs_grid/bk_line/use_bk_line.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4ec37fd8 on 2023-10-03 22:00:32 UTC
323aa12fe2 Jean*0001
0002 % input: krd, kfac, jprt
4ec37fd829 Jean*0003 % kfac = 0 : read-in 1 diagnostics output file containing both UVELMASS & VVELMASS
0004 % kfac > 0 : read-in 2 Time-ave output files ; kfac=1 with hFac in ; kfac=2 : add hFac here.
323aa12fe2 Jean*0005 krd=2; kfac=0; jprt=0;
0006
0007 if krd > 0,
4ec37fd829 Jean*0008 gDir='grid_files/'; rDir='output_files/'; nit=72000;
0009 %rDir='../res_s1y/'; nit=72360;
0010 %gDir=rDir;
0011 titexp=strrep(rDir(1:end-1),'_','\_'); tit0=['; it=',int2str(nit)];
0012 if kfac > 0,
0013 %-- time-ave output:
0014 uu=rdmds([rDir,'hUtave'],nit);
0015 vv=rdmds([rDir,'hVtave'],nit);
0016 if kfac==2,
0017 %- without included hFac:
0018 hFacW=rdmds([rDir,'hFacW']);
0019 hFacS=rdmds([rDir,'hFacS']);
0020 uu=hFacW.*uu; vv=hFacS.*vv;
0021 end
0022 else
0023 %-- diagnostics output:
0024 namFil='dynDiag';
0025 [v4d,dum,M]=rdmds([rDir,namFil],nit);
0026 eval(M);
0027 tmp=strcmp(fldList,'UVELMASS'); ju=find(tmp==1);
0028 tmp=strcmp(fldList,'VVELMASS'); jv=find(tmp==1);
0029 if isempty(ju) | isempty(jv),
0030 error(['At least 1 velocity Comp. missing from file: ',namFil]);
0031 end
0032 uu=squeeze(v4d(:,:,:,ju));
0033 vv=squeeze(v4d(:,:,:,jv));
0034 end
0035 n1h=size(uu,1); n2h=size(uu,2);
0036 if n1h == 6*n2h, nc=n2h;
0037 elseif n1h*6 == n2h, nc=n1h;
0038 else
0039 error([' var size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
323aa12fe2 Jean*0040 end
4ec37fd829 Jean*0041 nr=size(uu,3); nPg=nc*nc*6;
323aa12fe2 Jean*0042 end
0043
0044 if krd > 1,
4ec37fd829 Jean*0045 lDir='grid_files/';
28aaff1409 Jean*0046 if krd == 2,
0047 bk_lineF='isoLat_cs32_59';
0048 %- load broken lines (iso-lat) :
323aa12fe2 Jean*0049 % bkl_Ylat, bkl_Npts, bkl_Flg, bkl_IJuv, bkl_Xsg, bkl_Ysg, bkl_Zon
4ec37fd829 Jean*0050 load([lDir,bk_lineF]);
28aaff1409 Jean*0051 end
0052 if krd == 3,
0053 bk_lineF='bkl_AB_cs32';
0054 %- load broken line (great-circle arc AB):
0055 % ab_Npts, ab_Flg, ab_IJuv, ab_Xsg, ab_Ysg
4ec37fd829 Jean*0056 load([lDir,bk_lineF]);
0057 bkl_Ylat=0; bkl_Npts=ab_Npts;
28aaff1409 Jean*0058 bkl_Flg=ab_Flg'; bkl_IJuv=ab_IJuv'; bkl_Xsg=ab_Xsg'; bkl_Ysg=ab_Ysg';
0059 end
323aa12fe2 Jean*0060
0061 ufac=rem(bkl_Flg,2) ; vfac=fix(bkl_Flg/2) ;
0062
0063 %- load the grid dx,dy :
0c7ba451ba Jean*0064 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0065 % or ncdf=0 to load MDS (binary) grid-files :
0066 ncdf=0;
4ec37fd829 Jean*0067 G=load_grid(gDir,ncdf);
0c7ba451ba Jean*0068 dxg=G.dXg; dyg=G.dYg;
323aa12fe2 Jean*0069 fprintf([' load bk_line description from file=',bk_lineF, ...
0070 ' AND dxg,dyg files \n']);
0071 dxg=reshape(dxg,nPg,1); dyg=reshape(dyg,nPg,1);
0072 end
0073
0074 %- Atmospheric example:
0075 %delR=[10.E3, 25.E3, 30.E3, 20.E3, 15.E3]; % = delP in Pa
0076 %psiFac=1.e-9; g=9.81; % <- convert to 10^9 kg/s
0077 %- Oceanic example:
4ec37fd829 Jean*0078 %delR=[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]; % = delR in m
0079 delR=G.dRf;
0080 psiFac=-1.e-6; g=1; % <- convert to Sv. [10^6 m3/s]
323aa12fe2 Jean*0081
0082 %- compute the horizontal transport ut,vt :
0083 ut=reshape(uu,[nPg nr]); vt=reshape(vv,[nPg nr]);
0084 for k=nr:-1:1,
0085 ut(:,k)=psiFac*dyg.*ut(:,k);
0086 vt(:,k)=psiFac*dxg.*vt(:,k);
0087 end
4ec37fd829 Jean*0088 if n2h == 6*nc,
0089 %- to use this (older) bk-line file, long-vector is in "old-format" style, i.e., nc*6*nc:
0090 ut=reshape(ut,[nc nc 6 nr]); ut=permute(ut,[1 3 2 4]);
0091 vt=reshape(vt,[nc nc 6 nr]); vt=permute(vt,[1 3 2 4]);
0092 ut=reshape(ut,[nPg nr]); vt=reshape(vt,[nPg nr]);
0093 end
323aa12fe2 Jean*0094
0095 %- integrate along a broken-line (~ iso Latitude)
0096 ydim=length(bkl_Ylat); ylat=bkl_Ylat;
0097 vz=zeros(ydim,nr);
0098 for jl=1:ydim,
4ec37fd829 Jean*0099 if jl == jprt, fprintf(' jl= %2i , lat= %8.3f , Npts= %3i :\n', ...
323aa12fe2 Jean*0100 jl,ylat(jl),bkl_Npts(jl)); end
4ec37fd829 Jean*0101 if 1 < 0,
0102 %-- slower version (with debug print):
0103 for ii=1:bkl_Npts(jl),
0104 ij=bkl_IJuv(ii,jl);
0105 if jl == jprt,
0106 i = 1+rem(ij-1,6*nc); j=1+fix((ij-1)/nc/6); f=1+fix((i-1)/nc); i=1+rem(i-1,nc);
0107 fprintf(' no= %3i : Flg,I,J= %2i (%2i,%2i) %3i %3i (nf=%1i)\n', ...
0108 ii,bkl_Flg(ii,jl),ufac(ii,jl),vfac(ii,jl),i,j,f);
0109 end
0110 for k=1:nr,
0111 vz(jl,k)=vz(jl,k)+ufac(ii,jl)*ut(ij,k)+vfac(ii,jl)*vt(ij,k);
0112 end
323aa12fe2 Jean*0113 end
4ec37fd829 Jean*0114 else
0115 %-- faster version (no debug print):
0116 % Note: could reduce memory if permute k & jl loop and compute ut,vt within same k loop
0117 % and then can also cumulate stream-function within same single k loop
0118 ie=bkl_Npts(jl);
323aa12fe2 Jean*0119 for k=1:nr,
4ec37fd829 Jean*0120 vz(jl,k)=sum( ufac(1:ie,jl).*ut(bkl_IJuv(1:ie,jl),k) ...
0121 + vfac(1:ie,jl).*vt(bkl_IJuv(1:ie,jl),k) );
323aa12fe2 Jean*0122 end
0123 end
4ec37fd829 Jean*0124
323aa12fe2 Jean*0125 end
4ec37fd829 Jean*0126
323aa12fe2 Jean*0127 %- compute the meridional transport streamfunction:
0128 psi=zeros(ydim+2,nr+1);
0129 for j=1:ydim, for k=nr:-1:1,
0130 psi(j+1,k)=psi(j+1,k+1) + delR(k)*vz(j,k)/g ;
0131 end ; end
0132
0133 %-- make a graphe :
0134 yax=zeros(ydim+2,1);yax(2:ydim+1)=ylat;
28aaff1409 Jean*0135 if krd ==3, yax=[-10 0 10];
0136 else yax(1)=2*ylat(1)-ylat(2); yax(2+ydim)=2*ylat(ydim)-ylat(ydim-1); end
323aa12fe2 Jean*0137 zax=[0:nr];
0138 %---
4ec37fd829 Jean*0139 %ccp=[10:20:90]; cc=[-ccp(end:-1:1) 0 ccp];
0140 ccp=[5:5:20 30:10:90]; cc=[-ccp(end:-1:1) 0 ccp];
323aa12fe2 Jean*0141 [cs,h]=contour(yax,zax,psi',cc);clabel(cs); isoline0(h);
0142 if g==1, set(gca,'YDir','reverse'); end
0143 grid on;
4ec37fd829 Jean*0144 MxV=max(max(abs(psi)));
323aa12fe2 Jean*0145 Msurf=max(squeeze(abs(psi(:,1))));
0146 if g== 1,
4ec37fd829 Jean*0147 text(0,nr*1.07,sprintf('Max= %9.5g ; Surf-Max= %9.5g ', MxV, Msurf))
323aa12fe2 Jean*0148 title(['Merid. Transport [Sv] ; ',titexp,tit0]);
0149 else
0150 text(0,-0.5,sprintf('Max= %9.5g ; Surf-Max= %9.5g ', MxV, Msurf))
0151 title(['Merid. Transport [10^9 kg/s] ; ',titexp,tit0]);
0152 end
4ec37fd829 Jean*0153 %-----------------
323aa12fe2 Jean*0154
0155 return
