% ROMS_TRI_TIDES: Makes ROMS tidal forcing from triangle based tidal data % This script interpolates tidal elevation and current data from % triangle-based tidal data. This could be output from a Finite Element % Model such as ADCIRC or QUODDY or any data run through DELAUNAY triangulation. % % User must modify the file names, start time and tidal constituents % for the particular application in the user-modified section below. % % Rich Signell (rsignell@usgs.gov) % John Warner (jcwarner@usgs.gov) % % % Requires: TRI, T_TIDE, and TIDAL_ELLIPSE packages %%%%% BEGINNING of USER-MODIFIED SECTION %%%%% IPLOT=0; % 1 to make plots, 0 for no plots IWRITE=1; % 1 to write output to netcdf, 0 for no output % (1) Specify Grid File and Forcing File % (Create if needed) % Example: nc_create('your_file_name.nc',#x rho points,#y rho points,0,0) Gname='USeast_grd3.nc'; Fname='USeast_adc2001v2e_07-feb-2005.nc'; nc_create(Fname,896,336,0,0); % (2) Enter ROMS start time. This will be used to calculate the proper phase % for the tidal constituents. If you change the start time, this routine % must be run again to adjust the tidal phases. % YYYY MO DA HR MI SC g = [ 2005, 02, 07, 0, 0, 0]; disp(['Tidal Start Time =' datestr(g)]) % (3) Load TRI tidal constituent data load adcirc_ec2001v2e_fix.mat tri lon lat u v elev depth periods freq names % (4) Specify which tidal constituents to use. Look at "names" in the .mat file to % see what is available. Program will halt if these are not found. tides_to_use=['K1' 'O1' 'Q1' 'M2' 'S2' 'N2' 'K2']; %%%%% END of USER-MODIFIED SECTION %%%%% % make sure that tidal data is double tri=double(tri); lon=double(lon); lat=double(lat); u=double(u); v=double(v); elev=double(elev); depth=double(depth); [Ntide,ncol]=size(tides_to_use); for k=1:Ntide imatch=strmatch(tides_to_use(k,:),names); if isempty(imatch), disp(['Error: Did not find constituent ' tides_to_use(k,:) ' in .mat file']); return end indx(k)=imatch; end %--------------------------------------------------------------------- % Read in application grid data. %--------------------------------------------------------------------- rlon=nc_read(Gname,'lon_rho'); rlat=nc_read(Gname,'lat_rho'); [Lp,Mp]=size(rlon); rmask =nc_read(Gname,'mask_rho'); rangle=nc_read(Gname,'angle'); % do only water points % [iwater]=find(~rmask); % do all points iwater=[1:size(rmask,1)*size(rmask,2)]; % Calculate V,U,F arguments using Rich P's T_Tide. % These are the astronomical offsets % that allow you to go from amplitude and Greenwich phase to % real tidal predictions for a certain time period. See the % T_TIDE toolkit and the T_VUF.M program for more info. a=t_getconsts; %--------------------------------------------------------------------- % Interpolate tide data to application grid. %--------------------------------------------------------------------- % % Initialize arrays for interpolated tidal data. ei=ones(Lp,Mp,Ntide).*NaN; % elevation ui=ones(Lp,Mp,Ntide).*NaN; % u vi=ones(Lp,Mp,Ntide).*NaN; % v tmp=ones(Lp,Mp).*NaN; % Interpolate elevation and velocity from ADCIRC model data using % "triterp", linear-interpolation from triangulated grid. for kk=1:Ntide, k=indx(kk); disp(['Interpolating ' names(k,:)]); format loose tmp(iwater)=triterp(tri,lon,lat,elev(:,k),rlon(iwater),rlat(iwater)); ei(:,:,kk)=tmp; tmp(iwater)=triterp(tri,lon,lat,u(:,k),rlon(iwater),rlat(iwater)); ui(:,:,kk)=tmp; tmp(iwater)=triterp(tri,lon,lat,v(:,k),rlon(iwater),rlat(iwater)); vi(:,:,kk)=tmp; % find T_TIDE constituent names (a.name) that match % selected ADCIRC constituent names (names(tides_to_use)) iconst(kk)=strmatch(names(k,:), a.name); end clear tmp %*********************************************************************** % This is the call to t_vuf that % will correct the phase to be at the user specified time. Also, the amplitude % is corrected for nodal adjustment. start_year=g(1,1); start_month=g(1,2); start_day=g(1,3); start_hour=g(1,4); start_min=g(1,5); start_sec=g(1,6); omega=2*pi.*a.freq(iconst); %tidal frequencies in radians/hour % Reference latitude for 3rd order satellites (degrees) is % set to 55. You don't need to adjust this to your local latitude % It could also be set to NaN as in Xtide, with very little effect. % See T_VUF for more info. reflat=55; % T_VUF syntax for T_TIDE v1.1: % [vv,uu,ff]=t_vuf(datenum(start_year,start_month,start_day,start_hour,start_min,start_sec),iconst,reflat); % T_VUF syntax for T_TIDE v1.2: [vv,uu,ff]=t_vuf('nodal',datenum(start_year,start_month,start_day,start_hour,start_min,start_sec),iconst,reflat); %vv and uu are returned in cycles, so * by 360 to get degrees or * by 2 pi to get radians vv=vv*360; % convert vv to phase in degrees uu=uu*360; % convert uu to phase in degrees %adjust phase and ampl of tide to account for user specified time ei_phase=ones(Lp,Mp,Ntide).*NaN; % elevation ei_amplitude=ones(Lp,Mp,Ntide).*NaN; % elevation ui_phase=ones(Lp,Mp,Ntide).*NaN; % u vel ui_amplitude=ones(Lp,Mp,Ntide).*NaN; % u vel vi_phase=ones(Lp,Mp,Ntide).*NaN; % v vel vi_amplitude=ones(Lp,Mp,Ntide).*NaN; % v vel for k=1:Ntide; ei_phase(:,:,k) = angle(ei(:,:,k)).*180./pi - vv(k) - uu(k); % degrees ei_amplitude(:,:,k) = abs(ei(:,:,k)) .* ff(k); ui_phase(:,:,k) = angle(ui(:,:,k)).*180./pi - vv(k) - uu(k); % degrees ui_amplitude(:,:,k) = abs(ui(:,:,k)) .* ff(k); vi_phase(:,:,k) = angle(vi(:,:,k)).*180./pi - vv(k) - uu(k); % degrees vi_amplitude(:,:,k) = abs(vi(:,:,k)) .* ff(k); end Tide.period=periods(indx); Tide.component=a.name(iconst,:); Tide.Ephase = mod(ei_phase(:,:,:),360); Tide.Eamp = ei_amplitude(:,:,:); size(Tide.Eamp) %--------------------------------------------------------------------- % Convert tidal current amplitude and phase lag parameters to tidal % current ellipse parameters: Major axis, ellipticity, inclination, % and phase. Use "tidal_ellipse" package by Zhigang Xu %--------------------------------------------------------------------- ui_pha=ui_phase(:,:,:); ui_amp=ui_amplitude(:,:,:); vi_pha=vi_phase(:,:,:); vi_amp=vi_amplitude(:,:,:); [major,eccentricity,inclination,phase]=ap2ep(ui_amp,ui_pha,vi_amp,vi_pha); Tide.Cmax=major; Tide.Cmin=major.*eccentricity; Tide.Cangle=inclination; Tide.Cphase=phase; %clear ui ui_amp ui_pha vi vi_amp vi_pha ei %clear major eccentricity inclination phase %--------------------------------------------------------------------- % Plot results. %--------------------------------------------------------------------- if (IPLOT), for k = 1: Ntide figure subplot(1,2,1) pcolor(rlon,rlat,squeeze(Tide.Eamp(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Amplitude (m)'); title(strcat(Tide.component(k,:),' Tidal Component')); subplot(1,2,2) pcolor(rlon,rlat,squeeze(Tide.Ephase(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Phase Angle (degree)'); title(strcat(Tide.component(k,:),' Tidal Component')); end, for k=1:Ntide figure subplot(2,2,1) pcolor(rlon,rlat,squeeze(Tide.Cmax(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Major Axis amplitude (m/s)'); title(strcat(Tide.component(k,:),' Tidal Component')); subplot(2,2,2) pcolor(rlon,rlat,squeeze(Tide.Cmin(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Minor Axis amplitude (m/s)'); subplot(2,2,3) pcolor(rlon,rlat,squeeze(Tide.Cangle(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Inclination Angle (degrees)'); subplot(2,2,4) pcolor(rlon,rlat,squeeze(Tide.Cphase(:,:,k))); shading('interp'); colorbar; dasp(43); grid on; xlabel('Phase Angle(degrees)'); end, end, %--------------------------------------------------------------------- % Write out tide data into FORCING NetCDF file. %--------------------------------------------------------------------- if (IWRITE), [Vname,status]=wrt_tides(Fname,Tide); end