/[MITgcm]/MITgcm/verification/obcs_ctrl/input_ad/gendata.m
ViewVC logotype

Contents of /MITgcm/verification/obcs_ctrl/input_ad/gendata.m

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


Revision 1.1 - (show annotations) (download)
Wed Apr 20 19:18:26 2011 UTC (12 years, 11 months ago) by mmazloff
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, HEAD
Verification experiment using obcs controls

1 clear all
2 close all
3 ieee='b';accuracy='single';
4
5 dz = [500 500 500 500 500 500 500 500 ]
6 zf=-cumsum([0 dz]);
7 zc=(zf(1:end-1)+zf(2:end))/2;
8
9 Ho=-zf(length(zf));
10 nx=64;ny=64;nz = 8;
11
12 if 1
13 % set up input for model
14 % Flat bottom at z=-Ho
15 h=-Ho*ones(nx,ny);
16 %h(end,:)=0;%h(:,end)=0;%WALLS
17 fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
18
19 T = [20.,16.,12.,10., 9., 8., 7., 6.]
20 T2D = single(ones(64,1)*T);
21 % plot((T2D(1,:)),zc,'-ro','linewidth',2)
22 T3D = zeros(nx,ny,nz,'single');
23 for i = 1:64; T3D(i,:,:)= T2D; end
24
25 % T OBSERVED
26 fid = fopen('FinalThetaObs.bin','w',ieee);
27 fwrite(fid,T3D,accuracy);
28 fclose(fid)
29
30 %VPROFILE
31 x = linspace(-1,1,ny)';%x = [-1 linspace(-1,1,ny-2) 1]
32 V = single(1-(x.^2));
33 V2D =single(V*ones(1,nz))./10;
34 %V2D = zeros(64,nz,'single');
35 for i = 1:64
36 V3D(i,:,:) = V2D;
37 end
38
39 % Lets make a zonal jet solution
40 %INITIAL CONDITION
41 fid = fopen('Vini.bin','w',ieee);
42 fwrite(fid,V3D.*0,accuracy);
43 fclose(fid)
44 fid = fopen('Uini.bin','w',ieee);
45 fwrite(fid,V3D,accuracy);
46 fclose(fid)
47
48 %NORTHERN BOUNDARY
49 fid = fopen('Unbc.bin','w',ieee);
50 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
51 fclose(fid)
52 fid = fopen('Vnbc.bin','w',ieee);
53 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
54 fclose(fid)
55 fid = fopen('Snbc.bin','w',ieee);
56 fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
57 fclose(fid)
58 fid = fopen('Tnbc.bin','w',ieee);
59 fwrite(fid,T2D,accuracy);
60 fclose(fid)
61
62 %SOUTHERN BOUNDARY
63 fid = fopen('Usbc.bin','w',ieee);
64 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
65 fclose(fid)
66 fid = fopen('Vsbc.bin','w',ieee);
67 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
68 fclose(fid)
69 fid = fopen('Ssbc.bin','w',ieee);
70 fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
71 fclose(fid)
72 fid = fopen('Tsbc.bin','w',ieee);
73 fwrite(fid,T2D,accuracy);
74 fclose(fid)
75
76 %WESTERN BOUNDARY
77 fid = fopen('Uwbc.bin','w',ieee);
78 fwrite(fid,V2D,accuracy);
79 fclose(fid)
80 fid = fopen('Vwbc.bin','w',ieee);
81 fwrite(fid,zeros(ny,nz,accuracy),accuracy);
82 fclose(fid)
83 fid = fopen('Swbc.bin','w',ieee);
84 fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
85 fclose(fid)
86 fid = fopen('Twbc.bin','w',ieee);
87 fwrite(fid,T2D,accuracy);
88 fclose(fid)
89
90
91 %EASTERN BOUNDARY
92 fid = fopen('Uebc.bin','w',ieee);
93 fwrite(fid,V2D,accuracy);
94 fclose(fid)
95 fid = fopen('Vebc.bin','w',ieee);
96 fwrite(fid,zeros(ny,nz,accuracy),accuracy);
97 fclose(fid)
98 fid = fopen('Sebc.bin','w',ieee);
99 fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
100 fclose(fid)
101 fid = fopen('Tebc.bin','w',ieee);
102 fwrite(fid,T2D,accuracy);
103 fclose(fid)
104 % end of model setup if-then
105 end
106
107
108 if 0
109 % get vertical modes for obcs decomp
110 rhonil = 1035;g = 9.81;
111 RC = squeeze(rdmds('../GRID/RC'));
112 RF = squeeze(rdmds('../GRID/RF'));
113 DRF = squeeze(rdmds('../GRID/DRF'));
114 mask = squeeze(rdmds('../GRID/maskInC'));
115 RLOW=RF(end);
116 zmid = [0; (RC(1:end-1)+RC(2:end))/2;RLOW];
117 NZ = length(RC);
118
119 % all the dRhdz* come down to Nz
120 dRhodz_bar = rdmds('../run_fwd/dRhodz_5');
121 dRhodz_bar(mask==0)=0;dRhodz_bar(dRhodz_bar==0) = nan;
122 dRhodz_bar = squeeze(nanmean(nanmean(dRhodz_bar(32:33,32:33,:),1),2));
123 dRhodz_bar(1)=dRhodz_bar(2);dRhodz_bar(NZ+1)=dRhodz_bar(NZ);
124 Nz = (-g./rhonil.*dRhodz_bar).^0.5;
125
126 modesv = zeros(NZ,NZ,NZ);
127 % when only one depth, just have 1 for the mode
128 modesv(1,1,1) = 1.0;
129 % for more than 1 depth, solve eigenvalue problem
130 for k = 2:NZ;
131 % iz = vector of depth levels
132 iz = 1:k;
133 % print out a progress number
134 NZ-k,
135 % regularly spaced depths to the bottom of layer k;
136 % leave out the surface: VERT_FSFB2.m will supply it.
137 % 5 m spacing was selected as being small enough
138 zreg=-5:-5:RF(k+1);
139 Nzreg = interp1(zmid,Nz.^2,zreg,'linear',0);
140 [c2, F, G, N2, Pmid] = VERT_FSFB2(Nzreg,zreg);
141 %VERT_FSFB2.m adds a surface point
142 zreg = [0,zreg];
143 % sort from largest to smallest phase speed
144 [c2s indx] = sort(abs(c2),'descend');
145 % now interp back to our grid
146 for mds = 1:k
147 YI = interp1(zreg,F(:,indx(mds)),RC(iz),'linear',0);
148 modesv(iz,mds,k) = YI;
149 end
150 % have to weight by dz! Make a vector of fractional dz's
151 zwt = DRF(iz)/sum(DRF(iz),1);
152 %ensure first mode is barotropic (constant in depth)
153 avm1=sum(modesv(iz,1,k).*zwt,1);
154 modesv(iz,1,k)=avm1;
155 %make all modes orthogonal weighted by delta z
156 % use gram-schmidt leaving first one the same
157 for mds = 1:k-1
158 R = sum((modesv(iz,mds,k).^2).*zwt,1);
159 R2 = (modesv(iz,mds,k).*zwt)'*modesv(iz,mds+1:k,k)/R;
160 modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
161 modesv(iz,mds,k)*R2;
162 end
163 %All now orthognal, now normalize
164 for mds = 1:k
165 R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
166 if R < 1e-8
167 fprintf('Small R!! for mds = %2i and k = %2i\n',mds,k);
168 R = inf;
169 end
170 modesv(iz,mds,k) = modesv(iz,mds,k)./R;
171 end
172 end;% end of loop over depth level k
173 fid = fopen('baro_invmodes.bin','w','b');
174 fwrite(fid,modesv,'double');
175 fclose(fid)
176 if 1 %plot first 5 modes for deepest case
177 figure
178 clf;plot(modesv(:,[1:5],NZ),RC(:));
179 title('output modes at deepest point')
180 end
181 if 1 % test orthogonality
182 % do whole matrix (need to include dz!)
183 k = NZ;
184 cmm = (squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k]))'* ...
185 squeeze(modesv(iz,iz,k));
186 figure;imagesc(cmm);colorbar
187 title([num2str(k) ' mode orthogonality min, max diag: ' ...
188 num2str(min(diag(cmm))) ', ' ...
189 num2str(max(diag(cmm)))])
190 end
191 save baro_invmodes modesv
192 end
193

  ViewVC Help
Powered by ViewVC 1.1.22