/[MITgcm]/MITgcm_contrib/dgoldberg/depth_control_no_nsa/input/loadstate.m
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/depth_control_no_nsa/input/loadstate.m

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


Revision 1.1 - (show annotations) (download)
Thu Dec 7 23:28:44 2017 UTC (7 years, 7 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
*** empty log message ***

1 expdir = '../run';
2 iter=Inf;
3 g=mit_loadgrid;nx=g.nx; ny = g.ny;
4 grd = rdmnc(fullfile(expdir,'grid.*'),'Ro_surf','R_low','Depth');
5 g.depth=-grd.R_low;
6 g.ro_surf=-grd.Ro_surf;
7
8 DZ = repmat(reshape(g.dz,[1 1 g.nz]),[g.nx g.ny 1]);
9 s2d = rdmnc(fullfile(expdir,'surfDiag.*'),iter);
10 s3d = rdmnc(fullfile(expdir,'dynDiag.*'),iter);
11 e = s2d.ETAN(1:g.nx,1:g.ny,:,:);
12 t = s3d.THETA(1:g.nx,1:g.ny,:,:);
13 s = s3d.SALT(1:g.nx,1:g.ny,:,:);
14 uh = s3d.UVELMASS(1:g.nx,1:g.ny,:,:);
15 vh = s3d.VVELMASS(1:g.nx,1:g.ny,:,:);
16 % $$$ state = rdmnc(fullfile(expdir,'state.*'),iter);
17 % $$$ if isempty(state)
18 % $$$ [u,iter]=rdmds('U',iter);
19 % $$$ v=rdmds('V',iter);
20 % $$$ w=rdmds('W',iter);
21 % $$$ t=rdmds('T',iter);
22 % $$$ s=rdmds('S',iter);
23 % $$$ e=rdmds('Eta',iter);
24 % $$$ ph=rdmds('PH',iter);
25 % $$$ else
26 % $$$ u=state.U(1:g.nx,1:g.ny,:,:);
27 % $$$ v=state.V(1:g.nx,1:g.ny,:,:);
28 % $$$ w=state.W(1:g.nx,1:g.ny,:,:);
29 % $$$ t=state.Temp(1:g.nx,1:g.ny,:,:);
30 % $$$ s=state.S(1:g.nx,1:g.ny,:,:);
31 % $$$ e=state.Eta(1:g.nx,1:g.ny,:,:);
32 % $$$ phstate = rdmnc(fullfile(expdir,'phiHyd.*'),'phiHyd',iter);
33 % $$$ if size(e,3) > 1
34 % $$$ ph=phstate.phiHyd(1:g.nx,1:g.ny,:,[1 1:end]);
35 % $$$ else
36 % $$$ ph=phstate.phiHyd(1:g.nx,1:g.ny,:,:);
37 % $$$ end
38 % $$$ end
39 % $$$ ph0 = g.gravity*g.zc;
40 % $$$ msk = repmat(g.cmask,[1 1 1 size(u,4)]);
41 % $$$ %prs = g.rhonil*(ph+repmat(reshape(g.zc/g.gravity,[1 1 g.nz]),[nx ny 1 size(u,4)]));
42 % $$$ %prs = g.rhonil*(ph+repmat(reshape(e/g.gravity,[nx ny 1 size(u,4)]),[1 1 g.nz 1]));
43 % $$$ %prs = prs.*msk;
44 % $$$ hfw=change(g.hfacw,'==',NaN,0);
45 % $$$ hfs=change(g.hfacs,'==',NaN,0);
46 % $$$ uh = u;
47 % $$$ vh = v;
48 % $$$ for k=1:size(u,4)
49 % $$$ uh(:,:,:,k) = u(:,:,:,k).*hfw;
50 % $$$ vh(:,:,:,k) = v(:,:,:,k).*hfs;
51 % $$$ end
52 % $$$ [uc,vc,wc]=mit_velc(u,v,w);
53 [uhc,vhc]=mit_velc(uh,vh);
54
55 hfc=change(g.hfacc,'==',NaN,0);
56
57 xg=repmat(reshape(g.xg,[size(g.xg) 1]),[1 1 g.nz])/1000;
58 xc=repmat(reshape(g.xc,[size(g.xg) 1]),[1 1 g.nz])/1000;
59 yc=repmat(reshape(g.yc,[size(g.yg) 1]),[1 1 g.nz])/1000;
60 yg=repmat(reshape(g.yg,[size(g.yg) 1]),[1 1 g.nz])/1000;
61
62 fld=t;
63 %fld=uc;
64
65 % define section along main axis of PIG
66 [r,lo,la]=m_lldist([-102 -99.4],[-75 -75.3],20);
67 dd=m_lldist(lo,la);
68 xd = [0;cumsum(dd)];
69 for k=1:size(fld,3);
70 fldsec(:,k) = interp2(g.lonc,g.latc,(fld(:,:,k))',lo,la,'nearest');
71 end
72 depsec = -interp2(g.lonc,g.latc,g.depth',lo,la,'nearest');
73 rossec = -interp2(g.lonc,g.latc,g.ro_surf',lo,la,'nearest');
74
75 % bottom
76 landcolor = [1 1 1]*.7;
77 icecolor = [1 1 1]*1;
78
79 nn=length(xd);
80 ix = 1:nn; i=nn; k=0; while i > 1; ix=ix([1:i i i+1:nn+k]); i=i-1; k=k+1; end
81 xb = xd(ix);
82 zb = depsec([1 ix(1:end-1)]);
83 zb(end) = zb(1);
84 % top
85 zt = rossec([1 ix(1:end-1)]);
86 zt(end) = zt(1);
87
88 cax = [min(sq(fldsec(:))) max(sq(fldsec(:)))];
89 flev = cax(1):diff(cax)/20:cax(2);
90 h=pcol(xd,-g.zg,sq(fldsec)'); shading flat
91 hold on
92 fhb=fill([xb(1) xb' xb(end)],[-1000 zb -1000], ...
93 landcolor,'edgecolor','none');
94 fht=fill(xb,zt,icecolor,'edgecolor','k');
95 % $$$ fhb = plot(xb(:),sq(zb),'k');
96 % $$$ fht = plot(xb(:),sq(zt),'k');
97 hold off
98 caxis(cax);colorbar
99
100 pause;
101
102 m_proj('lambert','lon',[min(g.xg(:)) max(g.xg(:))+.125], ...
103 'lat',[min(g.yg(:)) max(g.yg(:))+.125*cos(max(g.yg(:))*pi/180)])
104 [xx,yy]=m_ll2xy(g.xg,g.yg);
105 [tfld,bfld] = topbot(g,fld);
106 clf
107 pcol(xx',yy',sq(tfld)');
108 m_grid;
109 hold on
110 sdz=sum(DZ.*hfc,3);
111 m_quiver(g.xc',g.yc',(sum(uhc.*DZ,3)./sdz)',(sum(vhc.*DZ,3)./sdz)',4,'k')
112 m_contour(g.lonc,g.latc,grd.Ro_surf',[-1 -1],'k')
113 hold off
114 %caxis([-2 -1.9])
115 colorbar
116
117 return
118 clf
119 pcol(xx',yy',sq(sq(grd.Ro_surf)-sq(grd.R_low))')
120 hold on
121 m_contour(g.xc',g.yc',(sq(grd.Ro_surf)-sq(grd.R_low))',[0:100:1000],'k')
122 caxis([0 700])
123 colorbar
124 m_grid
125
126 [x,y,z]=meshgrid(g.long,g.latg,g.zg);
127 sx = [-102:.2:-100];
128 sy = [-75.3:.1:-74.6];
129 sz = [-15 -500];
130 slice(x,y,-z,permute(sq(fld),[2 1 3]),sx,sy,sz);
131 shading flat
132 camlight
133 %caxis([-2 -1.88])
134
135 [ut,ub] = topbot(g,uhc./g.hfacc);
136 [vt,vb] = topbot(g,vhc./g.hfacc);
137 sdz=sum(DZ.*hfc,3);
138 uu=(sum(uhc.*DZ,3)./sdz);
139 vv=(sum(vhc.*DZ,3)./sdz);
140
141 pcol(xx',yy',-sq(s2d.SHIfwFlx)'/1000*31104000);
142 m_grid;
143 hold on
144 m_quiver(g.xc',g.yc',uu',vv',2,'k')
145 m_contour(g.lonc,g.latc,grd.Ro_surf',[-1 -1],'k')
146 hold off
147 %caxis([-2 -1.9])
148 colorbar

  ViewVC Help
Powered by ViewVC 1.1.22