Back to home page

MITgcm

 
 

    


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