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

Annotation 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 - (hide 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 dgoldberg 1.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