/[MITgcm]/MITgcm_contrib/high_res_cube/matlab/use_psiLine.m
ViewVC logotype

Contents of /MITgcm_contrib/high_res_cube/matlab/use_psiLine.m

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


Revision 1.1 - (show annotations) (download)
Wed Mar 28 18:07:27 2007 UTC (18 years, 3 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
Files contributed by JM for computation of barotropic streamfunction on cs510.

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

  ViewVC Help
Powered by ViewVC 1.1.22