%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Make plot from the results of the SHOREFACE test case % % Further Information: % http://www.crocoagrif.org/ % % This file is part of CROCOTOOLS % % CROCOTOOLS is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation; either version 2 of the License, % or (at your option) any later version. % % CROCOTOOLS is distributed in the hope that it will be useful, but % WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, % MA 02111-1307 USA % % Ref: Penven, P., L. Debreu, P. Marchesiello and J.C. McWilliams, % Application of the ROMS embedding procedure for the Central % California Upwelling System, Ocean Modelling, 2006. % % Patrick Marchesiello, IRD 2013 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all %================== User defined parameters =========================== % % --- model params --- % fname = 'shoreface_his.nc'; % croco file name yindex = 1; % y index makepdf = 0; % make pdf file % %====================================================================== if exist('subplot_tight')==2, mysubplot = str2func('subplot_tight'); else mysubplot = str2func('subplot'); end % --------------------------------------------------------------------- % --- get grid from numerical model --- % --------------------------------------------------------------------- nc=netcdf(fname,'r'); tindex=length(nc{'scrum_time'}(:)); % reads last record % % horizontal grid hr=squeeze(nc{'h'}(yindex,:)); xindex=1; hr=hr(xindex:end); L=length(hr); xr=squeeze(nc{'x_rho'}(yindex,xindex:end)); yr=squeeze(nc{'y_rho'}(yindex,xindex:end)); dx=xr(2)-xr(1); % % vertical grid N=length(nc('s_rho')); theta_s=nc.theta_s(:); theta_b=nc.theta_b(:); hc=nc.hc(:); zeta=squeeze(nc{'zeta'}(tindex,yindex,xindex:end)); zr=squeeze(zlevs(hr,zeta,theta_s,theta_b,hc,N,'r',2)); dzr=zr(2:end,:)-zr(1:end-1,:); % ---> zw(2:N,:) zru=0.5*(zr(:,1:end-1)+zr(:,2:end)); dzru=zru(2:end,:)-zru(1:end-1,:); % ---> zwu(2:N,:) zw=squeeze(zlevs(hr,zeta,theta_s,theta_b,hc,N,'w',2)); dzw=zw(2:end,:)-zw(1:end-1,:); % ---> zr zwu=0.5*(zw(:,1:end-1)+zw(:,2:end)); dzwu=zwu(2:end,:)-zwu(1:end-1,:); % ---> zru % xr2d=repmat(xr,[N 1]); xw2d=repmat(xr,[N+1 1]); D=hr+zeta; D2d=repmat(D,[N 1]); % --------------------------------------------------------------------- % --- read/compute numerical model fields (index 1) --- % -------------------------------------------------------------------- time=nc{'scrum_time'}(tindex)/86400; zeta1=zeta; % ... zonal velocity ... ---> xu,zru u1=squeeze(nc{'u'}(tindex,:,yindex,xindex:end)); % ... meridional velocity ... ---> xr,zr v1=squeeze(nc{'v'}(tindex,:,yindex,xindex:end)); % ... vertical velocity ... ---> xr,zw w1=squeeze(nc{'w'}(tindex,:,yindex,xindex:end)); % ... temperature ... ---> xr,zr t1=squeeze(nc{'temp'}(tindex,:,yindex,xindex:end)); % ... vertical viscosity/diffusivity Akv=squeeze(nc{'AKv'}(tindex,:,yindex,xindex:end)); Akt=squeeze(nc{'AKt'}(tindex,:,yindex,xindex:end)); % ... wave setup ... sup=squeeze(nc{'sup'}(tindex,yindex,xindex:end)); % ... zonal Stokes dritf ---> xu,zru ust=squeeze(nc{'ust'}(tindex,:,yindex,xindex:end)); % ... meridional Stokes drift ... ---> xr,zr vst=squeeze(nc{'vst'}(tindex,:,yindex,xindex:end)); % ... vertical Stokes drift ... ---> xr,zw wst=squeeze(nc{'wst'}(tindex,:,yindex,xindex:end)); % eddy viscosity due to depth-induced wave breaking Akb=squeeze(nc{'Akb'}(tindex,:,yindex,xindex:end)); % eddy diffusivity due to primary waves Akw=squeeze(nc{'Akw'}(tindex,:,yindex,xindex:end)); close(nc) %============================================================ % --- plot --- %============================================================= Dcrit=0.2; xr=1.e-3*xr-1; xr2d=1.e-3*xr2d-1; xw2d=1.e-3*xw2d-1; u1(:,L)=u1(:,L-1); zeta1(D