| 1 |
krd=2; kgr=1; kwr=0; Tprt=0; |
| 2 |
|
| 3 |
% $Header: /u/gcmpack/MITgcm/utils/matlab/cs_grid/bk_line/use_psiLine.m,v 1.1 2005/09/15 16:46:28 jmc Exp $ |
| 4 |
% $Name: $ |
| 5 |
|
| 6 |
if krd > 0, |
| 7 |
Rac='grid_files/'; |
| 8 |
drF=squeeze(rdmds([Rac,'DRF']))'; nr=length(drF); |
| 9 |
|
| 10 |
%- to make a plot with grph_CS: |
| 11 |
fprintf(' load xc,yc,xg,yg file ...'); |
| 12 |
xcs=rdmds([Rac,'XC']); |
| 13 |
ycs=rdmds([Rac,'YC']); |
| 14 |
xcg=rdmds([Rac,'XG']); |
| 15 |
ycg=rdmds([Rac,'YG']); |
| 16 |
fprintf(' done\n'); |
| 17 |
nc=size(xcs,2); nPg=nc*nc*6; |
| 18 |
%- add the 2 missing corners: |
| 19 |
xcg=reshape(xcg,[nPg 1]); ycg=reshape(ycg,[nPg 1]); |
| 20 |
xcg(nPg+1)=xcg(1); ycg(nPg+1)=ycg(1+2*nc); |
| 21 |
xcg(nPg+2)=xcg(1+3*nc); ycg(nPg+2)=ycg(1); |
| 22 |
|
| 23 |
psiLineF=[Rac,'psiLine_N2S_cs',int2str(nc)]; |
| 24 |
%psiLineF=['psiLine_N2S_cs',int2str(nc)]; |
| 25 |
|
| 26 |
fprintf([' load psiLine from file=',psiLineF,' ...']); |
| 27 |
%- load : 'psi_N','psi_C','psi_P','psiUV','psi_T' : |
| 28 |
load(psiLineF); |
| 29 |
fprintf(' done\n'); |
| 30 |
%- load the grid dx,dy , convert to 10^6 m : |
| 31 |
%dxg=rdmds([Rac,'DXG']); dxg=dxg*1.e-6; |
| 32 |
%dyg=rdmds([Rac,'DYG']); dyg=dyg*1.e-6; |
| 33 |
|
| 34 |
end |
| 35 |
|
| 36 |
if krd > 1, |
| 37 |
%-- load Vertically Integrated transport : |
| 38 |
rac='output_files/'; |
| 39 |
%nit=72020; |
| 40 |
%uu=rdmds([rac,'hUtave'],nit); |
| 41 |
%vv=rdmds([rac,'hVtave'],nit); |
| 42 |
namf=[rac,'UintTrsYr.bin']; |
| 43 |
fprintf([' read file: ',namf,' ...']); |
| 44 |
ut=rdda(namf,[nc*6*nc 1],1,'real*8','b'); |
| 45 |
fprintf(' done\n'); |
| 46 |
namf=[rac,'VintTrsYr.bin']; |
| 47 |
fprintf([' read file: ',namf,' ...']); |
| 48 |
vt=rdda(namf,[nc*6*nc 1],1,'real*8','b'); |
| 49 |
fprintf(' done\n'); |
| 50 |
end |
| 51 |
|
| 52 |
%- compute detph integrated flow : |
| 53 |
%delR=[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]; % = delR in m |
| 54 |
%ut=reshape(uu,nPg,nr); vt=reshape(vv,nPg,nr); |
| 55 |
%dxg=reshape(dxg,nPg,1); dyg=reshape(dyg,nPg,1); |
| 56 |
%var=dyg*delR; ut=var.*ut; ut=sum(ut,2); |
| 57 |
%var=dxg*delR; vt=var.*vt; vt=sum(vt,2); |
| 58 |
|
| 59 |
%- compute Barotropic Stream Function : |
| 60 |
psiNx=size(psi_C,1);psiNy=size(psi_C,2); nPc2=nPg+2; |
| 61 |
psiBa=zeros(nPc2,1); |
| 62 |
ufac=rem(psi_T,2) ; vfac=fix(psi_T/2) ; |
| 63 |
fprintf(' start computing stream-function ...'); |
| 64 |
if Tprt, TimeT0=clock; end |
| 65 |
for j=2:psiNy, |
| 66 |
for i=1:psi_N(j), |
| 67 |
i1=psi_C(i,j); i0=psi_P(i,j); i2=psiUV(i,j); |
| 68 |
psiBa(i1)=psiBa(i0)+ufac(i,j)*ut(i2)+vfac(i,j)*vt(i2); |
| 69 |
end |
| 70 |
end |
| 71 |
if Tprt, TimeT1=clock; end |
| 72 |
fprintf(' done\n'); |
| 73 |
if Tprt, fprintf(' Total time: %9.6f \n',etime(TimeT1,TimeT0) ); end |
| 74 |
|
| 75 |
if kwr==1, |
| 76 |
namfil='Psi_Horiz.bin'; |
| 77 |
fprintf([' write to file: ',namfil,' ...']); |
| 78 |
fid=fopen(namfil,'w','b'); fwrite(fid,psiBa,'real*8'); fclose(fid); |
| 79 |
fprintf(' done\n'); |
| 80 |
end |
| 81 |
|
| 82 |
if kgr > 0, |
| 83 |
shift=0; ccB=[0 0]; cbV=3; AxBx=[-180 180 -90 90]; kEnv=1; |
| 84 |
ccB=[-240 240]; |
| 85 |
%AxBx=[-210 210 30 90]; |
| 86 |
grph_CSz(psiBa,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv) |
| 87 |
end |