Warning, /utils/matlab/cs_grid/use_isoLat_bkl.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
4ec37fd829 Jean*0001
0002 % input: krd, kfac, jprt
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.
0005 krd=2; kfac=0; jprt=0;
0006
0007 if krd > 0,
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 !']);
0040 end
0041 nr=size(uu,3); nPg=nc*nc*6;
0042 end
0043
0044 if krd > 1,
0045 lDir='grid_files/';
0046 if krd == 2,
0047 bk_lineF='isoLat_59_cs32';
0048 %bk_lineF='isoLat_89_cs32';
0049 %- load broken lines (iso-lat) :
0050 % bkl_yLat, bkl_Npts, bkl_Flag, bkl_IJuv, bkl_tagC
0051 load([lDir,bk_lineF]);
0052 end
0053 if krd == 3,
0054 %-- new method has not been applied yet to great-circle broken line:
0055 error(' transport through great-circle line (krd=3) not supported here');
0056 bk_lineF='bkl_AB_cs32';
0057 %- load broken line (great-circle arc AB):
0058 % ab_Npts, ab_Flg, ab_IJuv, ab_Xsg, ab_Ysg
0059 load([lDir,bk_lineF]);
0060 bkl_yLat=0; bkl_Npts=ab_Npts;
0061 bkl_Flag=ab_Flg'; bkl_IJuv=ab_IJuv'; bkl_Xsg=ab_Xsg'; bkl_Ysg=ab_Ysg';
0062 end
0063
0064 ufac=rem(bkl_Flag,2) ; vfac=fix(bkl_Flag/2) ;
0065
0066 %- load the grid dx,dy :
0067 %- set ncdf=1 to load MNC (NetCDF) grid-files ;
0068 % or ncdf=0 to load MDS (binary) grid-files :
0069 ncdf=0;
0070 G=load_grid(gDir,ncdf);
0071 dxg=G.dXg; dyg=G.dYg;
0072 fprintf([' load bk_line description from file=',bk_lineF, ...
0073 ' AND dxg,dyg files \n']);
0074 dxg=reshape(dxg,nPg,1); dyg=reshape(dyg,nPg,1);
0075 end
0076
0077 %- Atmospheric example:
0078 %delR=[10.E3, 25.E3, 30.E3, 20.E3, 15.E3]; % = delP in Pa
0079 %psiFac=1.e-9; g=9.81; % <- convert to 10^9 kg/s
0080 %- Oceanic example:
0081 %delR=[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]; % = delR in m
0082 delR=G.dRf;
0083 psiFac=-1.e-6; g=1; % <- convert to Sv. [10^6 m3/s]
0084
0085 %- compute the horizontal transport ut,vt :
0086 ut=reshape(uu,[nPg nr]); vt=reshape(vv,[nPg nr]);
0087 for k=nr:-1:1,
0088 ut(:,k)=psiFac*dyg.*ut(:,k);
0089 vt(:,k)=psiFac*dxg.*vt(:,k);
0090 end
0091 if n1h == 6*nc,
0092 %- to use bk-line, long-vector needs to be "compact" format compatible:
0093 ut=reshape(ut,[nc 6 nc nr]); ut=permute(ut,[1 3 2 4]);
0094 vt=reshape(vt,[nc 6 nc nr]); vt=permute(vt,[1 3 2 4]);
0095 ut=reshape(ut,[nPg nr]); vt=reshape(vt,[nPg nr]);
0096 end
0097
0098 %- integrate along a broken-line (~ iso Latitude)
0099 ydim=length(bkl_yLat); ylat=bkl_yLat;
0100 vz=zeros(ydim,nr);
0101 for jl=1:ydim,
0102 if jl == jprt, fprintf(' jl= %2i , lat= %8.3f , Npts= %3i :\n', ...
0103 jl,ylat(jl),bkl_Npts(jl)); end
0104 if 1 < 0,
0105 %-- slower version:
0106 for ii=1:bkl_Npts(jl),
0107 ij=bkl_IJuv(ii,jl);
0108 for k=1:nr,
0109 vz(jl,k)=vz(jl,k)+ufac(ii,jl)*ut(ij,k)+vfac(ii,jl)*vt(ij,k);
0110 end
0111 end
0112 else
0113 %-- faster version:
0114 % Note: could reduce memory if permute k & jl loop and compute ut,vt within same k loop
0115 % and then can also cumulate stream-function within same single k loop
0116 ie=bkl_Npts(jl);
0117 for k=1:nr,
0118 vz(jl,k)=sum( ufac(1:ie,jl).*ut(bkl_IJuv(1:ie,jl),k) ...
0119 + vfac(1:ie,jl).*vt(bkl_IJuv(1:ie,jl),k) );
0120 end
0121 end
0122
0123 end
0124
0125 %- compute the meridional transport streamfunction:
0126 psi=zeros(ydim+2,nr+1);
0127 for j=1:ydim, for k=nr:-1:1,
0128 psi(j+1,k)=psi(j+1,k+1) + delR(k)*vz(j,k)/g ;
0129 end ; end
0130
0131 %-- make a graphe :
0132 yax=zeros(ydim+2,1);yax(2:ydim+1)=ylat;
0133 if krd ==3, yax=[-10 0 10];
0134 else yax(1)=2*ylat(1)-ylat(2); yax(2+ydim)=2*ylat(ydim)-ylat(ydim-1); end
0135 zax=[0:nr];
0136 %---
0137 %ccp=[10:20:90]; cc=[-ccp(end:-1:1) 0 ccp];
0138 ccp=[5:5:20 30:10:90]; cc=[-ccp(end:-1:1) 0 ccp];
0139 [cs,h]=contour(yax,zax,psi',cc);clabel(cs); isoline0(h);
0140 if g==1, set(gca,'YDir','reverse'); end
0141 grid on;
0142 MxV=max(max(abs(psi)));
0143 Msurf=max(squeeze(abs(psi(:,1))));
0144 if g== 1,
0145 text(0,nr*1.07,sprintf('Max= %9.5g ; Surf-Max= %9.5g ', MxV, Msurf))
0146 title(['Merid. Transport [Sv] ; ',titexp,tit0]);
0147 else
0148 text(0,-0.5,sprintf('Max= %9.5g ; Surf-Max= %9.5g ', MxV, Msurf))
0149 title(['Merid. Transport [10^9 kg/s] ; ',titexp,tit0]);
0150 end
0151 %-----------------
0152
0153 return
