% this code takes a list of u,v node numbers, and calculates the transport % between these sections, which can be averaged over a tidal cycle. % % also calculates T/S fluxes, useing scheme that the model does, % except (NB!) it is not downwind. % clf clear all orient landscape %profile on -detail builtin %connect the data files to be ploted filehist='gom_history.cdf'; filedir='/home/pringle/fvcom/fvcom15feb'; [jnk,basedir]=system(['basename ' filedir]); basedir=basedir(1:end-1) nchist=netcdf([filedir '/outmedm/' filehist],'nowrite'); t=nchist{'t1'}; s=nchist{'s1'}; u=nchist{'u'}; v=nchist{'v'}; el=nchist{'el'}; thour=nchist{'thour'}; %define the time and depth to get zplot=-200; nt=237; dnt=4; ntvec=1:1:(length(thour)-dnt+1); ntvec=1:1:length(thour); %ntvec=1:1:round(length(thour)/2); %ntvec=280:1:length(thour); %ntvec=(30*8):1:length(thour); %ntvec=25*8:28*8; %ntvec=1:272; %load mesh, contour depth load mesh subplot(2,2,1) load depthgrd hold on contour(depthgrd.xvec/1e3,depthgrd.yvec/1e3,... depthgrd.dpth,-[-270 -200 -100:20:0],'k-') hold off title('depth') set(gcf,'renderer','painters') %define the node number of the u,v points to calculate the transport between %a triangle xp=[1200 1200 1300 1200]*1e3+0*100e3; yp=[-100 0 -100 -100]*1e3-0*100e3; %the NEC %xp=[1.2219 1.2758]*1e6; %yp=[-86.5877 -37.8308]*1e3; %gom boundaries xp=[1.2813 1.2822 1.2503 0.9779 0.8850]*1e6; yp=[88.4808 -59.0878 -90.0589 -190.2598 -133.7829]*1e3; xp=[1.2813 1.2822 1.2503 0.9779 0.8850 1.2813]*1e6; yp=[88.4808 -59.0878 -90.0589 -190.2598 -133.7829 88.4808]*1e3; xp=[1.2813 1.2822 1.2503 0.9969 0.8850 0.9769 1.0696 1.2813]*1e6; yp=[88.4808 -59.0878 -90.0589 -197.7346 -133.7829 -95.0220 -45.8056 88.4808]*1e3; %Scotian Shelf coastal sec %xp=[1.4608 1.4772 1.3694 1.3405]*1e6; %yp=[2.3025 1.9695 1.2307 1.4374]*1e5; %Maine sec %xp=[1.0905 1.1009 0.9414 0.9256]*1e6; %yp=[1.9183 1.7264 0.9044 1.0709]*1e5; %CMO 70m isobath coastal sec %xp=[8.8906 8.7947 7.7706 7.7284]*1e5; %yp=[-1.3533 -2.6082 -2.4625 -1.6823]*1e5; %Maine sec %xp=[1.0905 1.1009 0.9414 0.9256 0.99045 1.0905]*1e6; %yp=[1.9183 1.7264 0.9044 1.0709 1.62453 1.9183]*1e5-500e3; %ax=[0.8272 1.3180 -0.0849 0.2274 ]*1e3; %xp=mesh.uvnode([1218 1216 1005 1218],1)'; %yp=mesh.uvnode([1218 1216 1005 1218],2)'; for n=1:length(xp) dist=sqrt((xp(n)-mesh.uvnode(:,1)).^2+(yp(n)-mesh.uvnode(:,2)).^2); [jnk,indx]=min(dist); nodevec(n)=indx; end nodevec %load existing path calculation, and if it is not the same as define by %nodevec, calculate a new path. if ~exist('shortpath.mat','file') shortpath.nodevec=[nan]; save shortpath shortpath end load shortpath if isequal(nodevec,shortpath.nodevec) dpath=shortpath.dpath; nends=shortpath.nends; else dpath=[]; nends=[]; for n=1:length(nodevec)-1 nends(n)=length(dpath)+1; temppath=shortest_path(mesh,nodevec(n+1),nodevec(n)); if (n~=(length(nodevec)-1)) dpath=[dpath temppath(1:end-1)]; else dpath=[dpath temppath(1:end)]; end end nends(end+1)=length(dpath); shortpath.dpath=dpath; shortpath.nodevec=nodevec; shortpath.nends=nends; save shortpath shortpath end % hold on % plot(mesh.uvnode(dpath,1)/1e3,mesh.uvnode(dpath,2)/1e3,'r-',... % 'linewidth',2); % hold off % transport is calculated from node to node in dpath. Note that this leaves % half a control volume out of the endpoint calcs. In the following figure, % "s" stands for the scalers, and "v" for the velocity points. % % s-------s % / \ / \ The big "V" mark the two ends of the path % / V-\-v-/-V \ Note that the flux through half the control % / \ / \ volume at either end is not included in the flux % s-------s------ s calculation. % % A positive flux is in the direction 90 degrees clockwise from the path, % as define by the direction from the first point in dpath to the last. % % In this case, if the leftmost "V" is the first one, a positive transport % would be down the page, and at right angles to the lines connecting the % "v" points. %calculate shared edges of adjoining triangles snodes_aft=zeros(length(dpath),2); snodes_for=zeros(length(dpath),2); for nnp=1:length(dpath) np=dpath(nnp); if (nnp~=1) snodes_aft(nnp,:)=... intersect(mesh.trinodes(np,:),mesh.trinodes(dpath(nnp-1),:)); end if (nnp~=length(dpath)) snodes_fore(nnp,:)=... intersect(mesh.trinodes(np,:),mesh.trinodes(dpath(nnp+1),:)); end end %plot sections on map colvec='rgbcym:' hold on for nsec=1:length(nends)-1 plot(mesh.uvnode(dpath(nends(nsec):nends(nsec+1)),1)/1e3,... mesh.uvnode(dpath(nends(nsec):nends(nsec+1)),2)/1e3,... [colvec(nsec)],... 'linewidth',3); end hold off %first, loop over all time steps overwhich the flux is to be calculated sflux=zeros(length(ntvec),length(dpath)); tflux=zeros(length(ntvec),length(dpath)); vflux=zeros(length(ntvec),length(dpath)); ntrip=0; for nt=ntvec nt ntrip=ntrip+1; %get all the data at that time step. tchunk=squeeze(t(nt,:,:)); schunk=squeeze(s(nt,:,:)); uchunk=squeeze(u(nt,:,:)); vchunk=squeeze(v(nt,:,:)); elchunk=squeeze(el(nt,:)); time=thour(nt); %loop over the path for nnp=1:length(dpath) np=dpath(nnp); %calculate depth averaged u,v, and calculate the free surface averaged %onto the location of the u,v point ua=mean(uchunk(:,np)); va=mean(vchunk(:,np)); uas=mean(uchunk(:,np).*(1/3).*(schunk(:,mesh.trinodes(np,1))+... schunk(:,mesh.trinodes(np,2))+schunk(:,mesh.trinodes(np,3)))); vas=mean(vchunk(:,np).*(1/3).*(schunk(:,mesh.trinodes(np,1))+... schunk(:,mesh.trinodes(np,2))+schunk(:,mesh.trinodes(np,3)))); uat=mean(uchunk(:,np).*(1/3).*(tchunk(:,mesh.trinodes(np,1))+... tchunk(:,mesh.trinodes(np,2))+tchunk(:,mesh.trinodes(np,3)))); vat=mean(vchunk(:,np).*(1/3).*(tchunk(:,mesh.trinodes(np,1))+... tchunk(:,mesh.trinodes(np,2))+tchunk(:,mesh.trinodes(np,3)))); el1=sum(elchunk(mesh.trinodes(np,:)))/3; h1=mesh.uvdepth(np); cent_x=mesh.uvnode(np,1); cent_y=mesh.uvnode(np,2); %except for the first node on the list, calculate transport across %line from the midpoint of a side of a t/s triangle node to the u,v point if (nnp~=1) %calculate midpoint of t/s triangle edge between last u,v node %and this one snodes=snodes_aft(nnp,:); half_x=0.5*(mesh.nodexy(snodes(1),1)+mesh.nodexy(snodes(2),1)); half_y=0.5*(mesh.nodexy(snodes(1),2)+mesh.nodexy(snodes(2),2)); if (nt==ntvec(1)) hold on plot([cent_x half_x]/1e3,[cent_y half_y]/1e3,'k') hold off end %calculate flux theta=atan2(cent_y-half_y,cent_x-half_x); ctheta=cos(theta); stheta=sin(theta); un=va*ctheta-ua*stheta; uns=vas*ctheta-uas*stheta; unt=vat*ctheta-uat*stheta; vflux(ntrip,nnp)=vflux(ntrip,nnp) + (el1+h1)*un*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); sflux(ntrip,nnp)=vflux(ntrip,nnp) + (el1+h1)*uns*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); tflux(ntrip,nnp)=tflux(ntrip,nnp) + (el1+h1)*unt*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); end %except for the last node on the list, calculate transport across %line from the u,v point to the midpoint of a side of a t/s triangle node if (nnp~=length(dpath)) %calculate midpoint of t/s triangle edge between last u,v node %and this one snodes=snodes_fore(nnp,:); half_x=0.5*(mesh.nodexy(snodes(1),1)+mesh.nodexy(snodes(2),1)); half_y=0.5*(mesh.nodexy(snodes(1),2)+mesh.nodexy(snodes(2),2)); if (nt==ntvec(1)) hold on plot([cent_x half_x]/1e3,[cent_y half_y]/1e3,'y') hold off end %calculate flux, note how angle calculation differs from above theta=atan2(half_y-cent_y,half_x-cent_x); ctheta=cos(theta); stheta=sin(theta); un=va*ctheta-ua*stheta; uns=vas*ctheta-uas*stheta; unt=vat*ctheta-uat*stheta; vflux(ntrip,nnp)=vflux(ntrip,nnp) + (el1+h1)*un*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); sflux(ntrip,nnp)=vflux(ntrip,nnp) + (el1+h1)*uns*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); tflux(ntrip,nnp)=tflux(ntrip,nnp) + (el1+h1)*unt*... sqrt((half_x-cent_x).^2+(half_y-cent_y).^2); end end end %filter time series? nfilt=8; nfvec=nfilt/2:(size(sflux,1)-nfilt/2); %this is done so that all of the points are included in a section nends(end)=nends(end)+1; subplot(2,2,2) for nsec=1:length(nends)-1; flux=real(boxcar_filt(sum(vflux(:,nends(nsec):(nends(nsec+1)-1)),2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),[ colvec(nsec)],... 'linewidth',2) hold on end flux=real(boxcar_filt(sum(vflux,2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),'-k','linewidth',2) hold off ylabel('volume flux') xlabel('tidal cycle') grid on subplot(2,2,3) for nsec=1:length(nends)-1; flux=real(boxcar_filt(sum(tflux(:,nends(nsec):(nends(nsec+1)-1)),2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),[ colvec(nsec)],... 'linewidth',2) hold on end flux=real(boxcar_filt(sum(tflux,2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),'-k','linewidth',2) hold off ylabel('temp flux, \int(u_nT)dxdz') xlabel('tidal cycle') grid on subplot(2,2,4) for nsec=1:length(nends)-1; flux=real(boxcar_filt(sum(sflux(:,nends(nsec):(nends(nsec+1)-1)),2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),[ colvec(nsec)],... 'linewidth',2) hold on end flux=real(boxcar_filt(sum(sflux,2),nfilt)); plot(thour(ntvec(nfvec))/12.42,flux(nfvec),'-k','linewidth',2) hold off ylabel('salt flux, \int(u_nS)dxdz') xlabel('tidal cycle') grid on suptitle(['on edge transport for ' basedir]) % profile report