%======================================================= % % $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/high_res_cube/eddy_flux/check_UV_STB.m,v 1.1 2004/08/16 21:10:22 edhill Exp $ % % Ed Hill % % Check the (u'b'),(v'b') values using the linear EoS: % % v'b' =approx= g * Alpha * v't' - g * Beta * v's' % % where: % g = 9.81 % Alpha = 2.10^{-4} % Beta = 7.10^{-4} clear all ; close all fnam = [ 'U_'; 'V_'; 'S_'; 'T_'; 'B_'; 'US'; 'UT'; 'UB'; 'VS'; 'VT'; 'VB' ]; vnam = lower(fnam); for i = 1:size(fnam,1) c = sprintf(['%sfid = ' ... 'fopen(''%s_tave_1996'',''r'',''ieee-be'');'], ... vnam(i,:),fnam(i,:) ); disp(c); eval(c); end delR = [ 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ... 10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ... 31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ... 93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ... 139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ... 341.50,364.50,387.50,410.50,433.50,456.50 ]; R = cumsum(delR) - 0.5*delR; nztot = 50; ne = 510; nps = 6 * ne * ne; iz = 1; for iz = 1:nztot offset = (iz - 1)*nps*4; for i = 1:size(fnam,1) c = sprintf('fseek(%sfid, offset, ''bof'');',vnam(i,:)); disp(c); eval(c); c = sprintf(['%s = reshape(fread(%sfid, nps,' ... '''real*4''),ne,ne,6);'],vnam(i,:),vnam(i,:)); disp(c); eval(c); end % Fill the (ne+1)*(ne+1) grid with T,S,B values ni = ne; nip1 = ni + 1; nj = ne; njp1 = nj + 1; tcubep1 = zeros( [ nip1 njp1 6 ] ); scubep1 = zeros( [ nip1 njp1 6 ] ); bcubep1 = zeros( [ nip1 njp1 6 ] ); for icf = 1:6 tcubep1(2:nip1,2:njp1,icf) = t_(:,:,icf); scubep1(2:nip1,2:njp1,icf) = s_(:,:,icf); bcubep1(2:nip1,2:njp1,icf) = b_(:,:,icf); end % Do the upwind-edge T exchanges tcubep1(1,2:njp1,1) = t_(ni:-1:1,nj,5); % - tcubep1(2:nip1,1,1) = t_(1:ni,nj,6); % - tcubep1(1,2:njp1,2) = t_(ni,1:nj,1); % - tcubep1(2:nip1,1,2) = t_(ni,nj:-1:1,6); % - tcubep1(1,2:njp1,3) = t_(ni:-1:1,nj,1); % - tcubep1(2:nip1,1,3) = t_(1:ni,nj,2); % - tcubep1(1,2:njp1,4) = t_(ni,1:nj,3); % - tcubep1(2:nip1,1,4) = t_(ni,nj:-1:1,2); % - tcubep1(1,2:njp1,5) = t_(ni:-1:1,nj,3); % - tcubep1(2:nip1,1,5) = t_(1:ni,nj,4); % - tcubep1(1,2:njp1,6) = t_(ni,1:nj,5); % - tcubep1(2:nip1,1,6) = t_(ni,nj:-1:1,4); % - % Do the upwind-edge S exchanges scubep1(1,2:njp1,1) = s_(ni:-1:1,nj,5); % - scubep1(2:nip1,1,1) = s_(1:ni,nj,6); % - scubep1(1,2:njp1,2) = s_(ni,1:nj,1); % - scubep1(2:nip1,1,2) = s_(ni,nj:-1:1,6); % - scubep1(1,2:njp1,3) = s_(ni:-1:1,nj,1); % - scubep1(2:nip1,1,3) = s_(1:ni,nj,2); % - scubep1(1,2:njp1,4) = s_(ni,1:nj,3); % - scubep1(2:nip1,1,4) = s_(ni,nj:-1:1,2); % - scubep1(1,2:njp1,5) = s_(ni:-1:1,nj,3); % - scubep1(2:nip1,1,5) = s_(1:ni,nj,4); % - scubep1(1,2:njp1,6) = s_(ni,1:nj,5); % - scubep1(2:nip1,1,6) = s_(ni,nj:-1:1,4); % - % Do the upwind-edge B exchanges bcubep1(1,2:njp1,1) = b_(ni:-1:1,nj,5); % - bcubep1(2:nip1,1,1) = b_(1:ni,nj,6); % - bcubep1(1,2:njp1,2) = b_(ni,1:nj,1); % - bcubep1(2:nip1,1,2) = b_(ni,nj:-1:1,6); % - bcubep1(1,2:njp1,3) = b_(ni:-1:1,nj,1); % - bcubep1(2:nip1,1,3) = b_(1:ni,nj,2); % - bcubep1(1,2:njp1,4) = b_(ni,1:nj,3); % - bcubep1(2:nip1,1,4) = b_(ni,nj:-1:1,2); % - bcubep1(1,2:njp1,5) = b_(ni:-1:1,nj,3); % - bcubep1(2:nip1,1,5) = b_(1:ni,nj,4); % - bcubep1(1,2:njp1,6) = b_(ni,1:nj,5); % - bcubep1(2:nip1,1,6) = b_(ni,nj:-1:1,4); % - % Get T values on the U,V grid points masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2)); maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1)); diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2); diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1); tonu = t_(:,:,:) + masku.*diffu; tonv = t_(:,:,:) + maskv.*diffv; % Get S values on the U,V grid points diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2); diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1); sonu = s_(:,:,:) + masku.*diffu; sonv = s_(:,:,:) + maskv.*diffv; % Get B values on the U grid points diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2); diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1); bonu = b_(:,:,:) + masku.*diffu; bonv = b_(:,:,:) + maskv.*diffv; % Calculate the primes uptp = ut - u_.*tonu; vptp = vt - v_.*tonv; upsp = us - u_.*sonu; vpsp = vs - v_.*sonv; upbp = ub - u_.*bonu; vpbp = vb - v_.*bonv; % Calculate the linear-EoS approximation Alpha = 2.0*10^-4; Beta = 7.0*10^-4; lin_upbp = 9.81.*( Alpha.*uptp - Beta.*upsp ); lin_vpbp = 9.81.*( Alpha.*vptp - Beta.*vpsp ); lns_upbp = 9.81.*( Alpha.*uptp ); ljs_upbp = 9.81.*( -Beta.*upsp ); % Quick plots figure(1), surf( vptp(:,:,6)), view(2), shading interp, colorbar, ... title('v''T''') figure(1), surf( upbp(:,:,1)), view(2), shading interp, colorbar, ... title('u''B'' from JMD95 EoS') figure(2), surf(lin_upbp(:,:,1)), view(2), shading interp, colorbar, ... title('u''B'' from Linear EoS') figure(3), surf(lns_upbp(:,:,1)), view(2), shading interp, colorbar, ... title('u''B'' by forgetting salt') figure(4), surf(ljs_upbp(:,:,1)), view(2), shading interp, colorbar, ... title('u''B'' with just salt') figure(3), surf( uptp(:,:,1)), view(2), shading interp, colorbar, ... title('u''T''') figure(4), surf( upsp(:,:,1)), view(2), shading interp, colorbar, ... title('u''S''') figure(3), surf( s_(:,:,1)), view(2), shading interp, colorbar, ... title('S') figure(4), surf( t_(:,:,1)), view(2), shading interp, colorbar, ... title('T') figure(3), surf( sonu(:,:,1)), view(2), shading interp, colorbar, ... title('S (on U points)') figure(4), surf( tonu(:,:,1)), view(2), shading interp, colorbar, ... title('T (on U points)') figure(3), surf( ut(:,:,1)), view(2), shading interp, colorbar, ... title('UT') figure(4), surf( us(:,:,1)), view(2), shading interp, colorbar, ... title('US') % Check bar(B) with bar(T),bar(S) p = 0.09998 * 9.81 * R(iz); for icf = 1:6 bcheck(:,:,icf) = ... eh3densjmd95( squeeze(s_(:,:,icf)), ... squeeze(t_(:,:,icf)), p ) * -9.81 / 1000; end figure(3), surf( b_(:,:,1)), view(2), shading interp, ... caxis([-10.1 -9.8]), colorbar, title('bar(B) from model') figure(4), surf( bcheck(:,:,1)), view(2), shading interp, ... caxis([-10.1 -9.8]), colorbar, title('bar(B) from bar(T),bar(S)') end