Warning, /verification/obcs_ctrl/input_ad/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit bce9e67d on 2011-04-20 19:18:26 UTC
bce9e67dfb Matt*0001 clear all
0002 close all
0003 ieee='b';accuracy='single';
0004
0005 dz = [500 500 500 500 500 500 500 500 ]
0006 zf=-cumsum([0 dz]);
0007 zc=(zf(1:end-1)+zf(2:end))/2;
0008
0009 Ho=-zf(length(zf));
0010 nx=64;ny=64;nz = 8;
0011
0012 if 1
0013 % set up input for model
0014 % Flat bottom at z=-Ho
0015 h=-Ho*ones(nx,ny);
0016 %h(end,:)=0;%h(:,end)=0;%WALLS
0017 fid=fopen('topog.box','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
0018
0019 T = [20.,16.,12.,10., 9., 8., 7., 6.]
0020 T2D = single(ones(64,1)*T);
0021 % plot((T2D(1,:)),zc,'-ro','linewidth',2)
0022 T3D = zeros(nx,ny,nz,'single');
0023 for i = 1:64; T3D(i,:,:)= T2D; end
0024
0025 % T OBSERVED
0026 fid = fopen('FinalThetaObs.bin','w',ieee);
0027 fwrite(fid,T3D,accuracy);
0028 fclose(fid)
0029
0030 %VPROFILE
0031 x = linspace(-1,1,ny)';%x = [-1 linspace(-1,1,ny-2) 1]
0032 V = single(1-(x.^2));
0033 V2D =single(V*ones(1,nz))./10;
0034 %V2D = zeros(64,nz,'single');
0035 for i = 1:64
0036 V3D(i,:,:) = V2D;
0037 end
0038
0039 % Lets make a zonal jet solution
0040 %INITIAL CONDITION
0041 fid = fopen('Vini.bin','w',ieee);
0042 fwrite(fid,V3D.*0,accuracy);
0043 fclose(fid)
0044 fid = fopen('Uini.bin','w',ieee);
0045 fwrite(fid,V3D,accuracy);
0046 fclose(fid)
0047
0048 %NORTHERN BOUNDARY
0049 fid = fopen('Unbc.bin','w',ieee);
0050 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
0051 fclose(fid)
0052 fid = fopen('Vnbc.bin','w',ieee);
0053 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
0054 fclose(fid)
0055 fid = fopen('Snbc.bin','w',ieee);
0056 fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
0057 fclose(fid)
0058 fid = fopen('Tnbc.bin','w',ieee);
0059 fwrite(fid,T2D,accuracy);
0060 fclose(fid)
0061
0062 %SOUTHERN BOUNDARY
0063 fid = fopen('Usbc.bin','w',ieee);
0064 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
0065 fclose(fid)
0066 fid = fopen('Vsbc.bin','w',ieee);
0067 fwrite(fid,zeros(nx,nz,accuracy),accuracy);
0068 fclose(fid)
0069 fid = fopen('Ssbc.bin','w',ieee);
0070 fwrite(fid,35.*ones(nx,nz,accuracy),accuracy);
0071 fclose(fid)
0072 fid = fopen('Tsbc.bin','w',ieee);
0073 fwrite(fid,T2D,accuracy);
0074 fclose(fid)
0075
0076 %WESTERN BOUNDARY
0077 fid = fopen('Uwbc.bin','w',ieee);
0078 fwrite(fid,V2D,accuracy);
0079 fclose(fid)
0080 fid = fopen('Vwbc.bin','w',ieee);
0081 fwrite(fid,zeros(ny,nz,accuracy),accuracy);
0082 fclose(fid)
0083 fid = fopen('Swbc.bin','w',ieee);
0084 fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
0085 fclose(fid)
0086 fid = fopen('Twbc.bin','w',ieee);
0087 fwrite(fid,T2D,accuracy);
0088 fclose(fid)
0089
0090
0091 %EASTERN BOUNDARY
0092 fid = fopen('Uebc.bin','w',ieee);
0093 fwrite(fid,V2D,accuracy);
0094 fclose(fid)
0095 fid = fopen('Vebc.bin','w',ieee);
0096 fwrite(fid,zeros(ny,nz,accuracy),accuracy);
0097 fclose(fid)
0098 fid = fopen('Sebc.bin','w',ieee);
0099 fwrite(fid,35.*ones(ny,nz,accuracy),accuracy);
0100 fclose(fid)
0101 fid = fopen('Tebc.bin','w',ieee);
0102 fwrite(fid,T2D,accuracy);
0103 fclose(fid)
0104 % end of model setup if-then
0105 end
0106
0107
0108 if 0
0109 % get vertical modes for obcs decomp
0110 rhonil = 1035;g = 9.81;
0111 RC = squeeze(rdmds('../GRID/RC'));
0112 RF = squeeze(rdmds('../GRID/RF'));
0113 DRF = squeeze(rdmds('../GRID/DRF'));
0114 mask = squeeze(rdmds('../GRID/maskInC'));
0115 RLOW=RF(end);
0116 zmid = [0; (RC(1:end-1)+RC(2:end))/2;RLOW];
0117 NZ = length(RC);
0118
0119 % all the dRhdz* come down to Nz
0120 dRhodz_bar = rdmds('../run_fwd/dRhodz_5');
0121 dRhodz_bar(mask==0)=0;dRhodz_bar(dRhodz_bar==0) = nan;
0122 dRhodz_bar = squeeze(nanmean(nanmean(dRhodz_bar(32:33,32:33,:),1),2));
0123 dRhodz_bar(1)=dRhodz_bar(2);dRhodz_bar(NZ+1)=dRhodz_bar(NZ);
0124 Nz = (-g./rhonil.*dRhodz_bar).^0.5;
0125
0126 modesv = zeros(NZ,NZ,NZ);
0127 % when only one depth, just have 1 for the mode
0128 modesv(1,1,1) = 1.0;
0129 % for more than 1 depth, solve eigenvalue problem
0130 for k = 2:NZ;
0131 % iz = vector of depth levels
0132 iz = 1:k;
0133 % print out a progress number
0134 NZ-k,
0135 % regularly spaced depths to the bottom of layer k;
0136 % leave out the surface: VERT_FSFB2.m will supply it.
0137 % 5 m spacing was selected as being small enough
0138 zreg=-5:-5:RF(k+1);
0139 Nzreg = interp1(zmid,Nz.^2,zreg,'linear',0);
0140 [c2, F, G, N2, Pmid] = VERT_FSFB2(Nzreg,zreg);
0141 %VERT_FSFB2.m adds a surface point
0142 zreg = [0,zreg];
0143 % sort from largest to smallest phase speed
0144 [c2s indx] = sort(abs(c2),'descend');
0145 % now interp back to our grid
0146 for mds = 1:k
0147 YI = interp1(zreg,F(:,indx(mds)),RC(iz),'linear',0);
0148 modesv(iz,mds,k) = YI;
0149 end
0150 % have to weight by dz! Make a vector of fractional dz's
0151 zwt = DRF(iz)/sum(DRF(iz),1);
0152 %ensure first mode is barotropic (constant in depth)
0153 avm1=sum(modesv(iz,1,k).*zwt,1);
0154 modesv(iz,1,k)=avm1;
0155 %make all modes orthogonal weighted by delta z
0156 % use gram-schmidt leaving first one the same
0157 for mds = 1:k-1
0158 R = sum((modesv(iz,mds,k).^2).*zwt,1);
0159 R2 = (modesv(iz,mds,k).*zwt)'*modesv(iz,mds+1:k,k)/R;
0160 modesv(iz,mds+1:k,k) = modesv(iz,mds+1:k,k) - ...
0161 modesv(iz,mds,k)*R2;
0162 end
0163 %All now orthognal, now normalize
0164 for mds = 1:k
0165 R = sqrt(sum((modesv(iz,mds,k).^2).*zwt,1));
0166 if R < 1e-8
0167 fprintf('Small R!! for mds = %2i and k = %2i\n',mds,k);
0168 R = inf;
0169 end
0170 modesv(iz,mds,k) = modesv(iz,mds,k)./R;
0171 end
0172 end;% end of loop over depth level k
0173 fid = fopen('baro_invmodes.bin','w','b');
0174 fwrite(fid,modesv,'double');
0175 fclose(fid)
0176 if 1 %plot first 5 modes for deepest case
0177 figure
0178 clf;plot(modesv(:,[1:5],NZ),RC(:));
0179 title('output modes at deepest point')
0180 end
0181 if 1 % test orthogonality
0182 % do whole matrix (need to include dz!)
0183 k = NZ;
0184 cmm = (squeeze(modesv(iz,iz,k)).*repmat(zwt,[1 k]))'* ...
0185 squeeze(modesv(iz,iz,k));
0186 figure;imagesc(cmm);colorbar
0187 title([num2str(k) ' mode orthogonality min, max diag: ' ...
0188 num2str(min(diag(cmm))) ', ' ...
0189 num2str(max(diag(cmm)))])
0190 end
0191 save baro_invmodes modesv
0192 end
0193
