clear %extractTPXOtides.m version 1.0. %Prepares a tidal forcing file for ROMS from TPXO tidal model (OSU). %Requires SNCtools and TPXO files h_tpxo72.nc, u_tpxo72.nc, and grid_tpxo71.nc %Plots require m_map toolbox and t_tide. %All those things are freely available on the web. %Phases are computed for time "t0", the instant that your simulation starts. %In this program t0 is a Matlab datenumber. The equivalent datenumber %used in ROMS is tR=t0-datenum(1968,5,23) (using Matlab function datenum). %If TIME_REF=-2 in sarom.in, TIDE_START=tR in sarom.in, and time (e.g. ocean_time)=tR in %the initial file, your tidal phase will be correct. I recommend also %that you set DSTART (in sarom.in) equal to tR so that the time stamp %in the standard output will be correct. %I'm reasonably confident that it's working correctly but let me know if %you find something wrong. Also, I'm not responsible for grounded oil %tankers, spoiled weddings on the beach, or anything else bad resulting from %the use of this program. On the other hand, if you like it, why not perform %a random act of kindness or something. Please contact me via the ROMS %forum if you find a bug or if I have left out a script from the zip file. %I have added 2D velocity output in an effort to encourage the ROMS %programmers to eventually eliminate the use of ellipses but they %(tide_Uamp, tide_Uphase, tide_Vamp, tide_Vphase) can be removed from the %script to neaten things up. For now ROMS doesn't use them anyway. %Note for rotated grids: velocities and ellipses are rotated in ROMS so they %are not rotated here. %J. Luick December 2011 via the ROMS Forum. %Modify the section between the asterisks for your given application. %************************************************************************************* TPXOfile_elev='h_tpxo72.nc'; %TPXO model file 1 %Must be on your path. TPXOfile_vel='u_tpxo72.nc'; %TPXO model file 1 %Download all three. TPXOfile_grid='grid_tpxo71.nc'; %TPXO model file 1 %from the TPXO website. workingdir='D:\MODELS\myroms\Projects\SAROM\DATA\'; %Windows version %workingdir='storage0/MODELS/myroms/Projects/SAROM/DATA/'; %Linux version ROMStitle='TPXO data for ROMS Model Run Dk'; ROMS_grid_file=[workingdir 'tidaltest_D_grid.nc']; ROMS_frc_file=[workingdir 'tidaltest_D_02JUL05.nc']; t0=datenum('02-Jul-2005 00:00'); %Time for which phase is computed. %Could be any HH:MM not necessarily 00:00 lengthSim=100; %Approximate anticipated length of model run (days) (for f & u) %ROMSnames={'K1' 'M2'}; %(for example) to use K1 and M2 only ROMSnames={'MM' 'MF' 'Q1' 'O1' 'P1' 'K1' 'N2' 'M2' 'S2' 'K2' 'MN4' 'M4' 'MS4'}; %to use all %Optional plotting stuff: sealevel_obs_file=[workingdir 'tidaltest_D_tideguage.dat']; %observations file (optional) %observations file to check timeseries. Use empty brackets if none available. %Column 1: datenumber %Column2: sealevel XtidePORT='Thevenard'; %Port to plot check timeseries (requires t_tide which uses Xtide) %Use empty quotes ('') to not plot an Xtide port. %Of course not all ports are in Xtide. LTZ=9.5; %Local standard time zone (+east of Greenwich; -west of Greenwich) (for Xtide plot only) %available choices (in TPXO): 'MM' 'MF' 'Q1' 'O1' 'P1' 'K1' 'N2' 'M2' 'S2' 'K2' 'MN4' 'M4' 'MS4' plot_timeseries_Lon=133.65; %Lon and Lat of ROMS rho-point in output file to plot plot_timeseries_Lat=-33; %a test timeseries MapPlotName='M2'; %Harmonic to plot in test timeseries ellipsePlotScale=1; %e.g. double to double the size of ellipses ellipsePlotStride=10; %e.g. make=10 to plot only each 10th ellipse %************************************************************************** %ROMS grid info lonR=nc_varget(ROMS_grid_file,'lon_rho'); latR=nc_varget(ROMS_grid_file,'lat_rho'); maskR=nc_varget(ROMS_grid_file,'mask_rho'); maskP=nc_varget(ROMS_grid_file,'mask_psi'); [M,L]=size(maskP); lonU=nc_varget(ROMS_grid_file,'lon_u'); latU=nc_varget(ROMS_grid_file,'lat_u'); lonV=nc_varget(ROMS_grid_file,'lon_v'); latV=nc_varget(ROMS_grid_file,'lat_v'); h=nc_varget(ROMS_grid_file,'h'); %Magic numbers for TPXO (Name, Solar Doodson Numbers, Ref. Phase, Speed degrees/hour) TPXOharmonics={ 'MM' '0' '1' '0' '-1' '0' '0' '0' '0.5443747' 'MF' '0' '2' '0' '0' '0' '0' '0' '1.0980331' 'Q1' '1' '-3' '1' '1' '0' '0' '270' '13.3986607' 'O1' '1' '-2' '1' '0' '0' '0' '270' '13.9430351' 'P1' '1' '0' '-1' '0' '0' '0' '270' '14.9589310' 'K1' '1' '0' '1' '0' '0' '0' '90' '15.0410690' 'N2' '2' '-3' '2' '1' '0' '0' '0' '28.4397297' 'M2' '2' '-2' '2' '0' '0' '0' '0' '28.9841042' 'S2' '2' '0' '0' '0' '0' '0' '0' '30.0000000' 'K2' '2' '0' '2' '0' '0' '0' '0' '30.0821381' 'MN4' '4' '-5' '4' '1' '0' '0' '0' '57.4238319' 'M4' '4' '-4' '4' '0' '0' '0' '0' '57.9682083' 'MS4' '4' '-2' '2' '0' '0' '0' '0' '58.9841042'}; %Magic numbers for ROMS Nharmonics=length(ROMSnames); ROMSharmonics=cell(Nharmonics,9); for n=1:Nharmonics m=find(strcmp(ROMSnames(n),TPXOharmonics(:,1))); ROMSharmonics(n,1:9)=TPXOharmonics(m,1:9); end ROMSperiods=360./str2double(ROMSharmonics(:,9)); % V,u,f (reference phase and nodal corrections) [fFac,uFac]=TPXOnodalfactors(t0+lengthSim/2,ROMSnames); Vdeg=Vphase(t0,ROMSharmonics); % Extract tide info from TPXO and put on rho grid zamp=zeros(Nharmonics,M+1,L+1); zpha=zeros(Nharmonics,M+1,L+1); uamp=zeros(Nharmonics,M+1,L+1); upha=zeros(Nharmonics,M+1,L+1); vamp=zeros(Nharmonics,M+1,L+1); vpha=zeros(Nharmonics,M+1,L+1); major=zeros(Nharmonics,M+1,L+1); eccentricity=zeros(Nharmonics,M+1,L+1); inclination=zeros(Nharmonics,M+1,L+1); phase=zeros(Nharmonics,M+1,L+1); for k=1:Nharmonics disp(['Extracting ',ROMSnames(k),' amplitudes']) zr=interpTPXO('hRe',ROMSnames(k),TPXOfile_elev,TPXOfile_grid,lonR,latR,maskR); zi=interpTPXO('hIm',ROMSnames(k),TPXOfile_elev,TPXOfile_grid,lonR,latR,maskR); ei=complex(zr,zi); zamp(k,:,:)=abs(ei).*fFac(k); zpha(k,:,:)=mod(-angle(ei)*180/pi-uFac(k)-Vdeg(k),360); disp(['Extracting ',ROMSnames(k),' u components']) ur=interpTPXO('URe',ROMSnames(k),TPXOfile_vel,TPXOfile_grid,lonR,latR,maskR); ui=interpTPXO('UIm',ROMSnames(k),TPXOfile_vel,TPXOfile_grid,lonR,latR,maskR); ei=complex(ur,ui); uamp(k,:,:)=abs(ei).*fFac(k); upha(k,:,:) =mod(-angle(ei)*180/pi-uFac(k)-Vdeg(k),360); disp(['Extracting ',ROMSnames(k),' v components']) vr=interpTPXO('VRe',ROMSnames(k),TPXOfile_vel,TPXOfile_grid,lonR,latR,maskR); vi=interpTPXO('VIm',ROMSnames(k),TPXOfile_vel,TPXOfile_grid,lonR,latR,maskR); ei=complex(vr,vi); vamp(k,:,:)=abs(ei).*fFac(k); vpha(k,:,:)=mod(-angle(ei)*180/pi-uFac(k)-Vdeg(k),360); [maj,ecc,inc,pha]=ap2ep(squeeze(uamp(k,:,:)),squeeze(upha(k,:,:)),... squeeze(vamp(k,:,:)),squeeze(vpha(k,:,:))); ecc(isnan(ecc))=0; %zero current results in e=NaN. ROMS crashes. (ps circle has e=0) major(k,:,:)=maj; eccentricity(k,:,:)=ecc; inclination(k,:,:)=inc; phase(k,:,:)=pha; end minor=major.*eccentricity; %Set up output netcdf forcing file nc_create_empty(ROMS_frc_file); %global atts %String version of ROMSnames (for global attribute only) c1=char(ROMSnames); [nr,nc]=size(c1); for i=1:Nharmonics ROMSnames_att((nc+1)*i-nc:(nc+1)*i)=[c1(i,1:nc) ' ']; end today=datestr(date,'yyyymmdd'); nc_attput(ROMS_frc_file,nc_global,'title',ROMStitle) nc_attput(ROMS_frc_file,nc_global,'Creation',today) nc_attput(ROMS_frc_file,nc_global,'grd_file',ROMS_grid_file) nc_attput(ROMS_frc_file,nc_global,'type','ROMS forcing file from TPXO') nc_attput(ROMS_frc_file,nc_global,'ini date (Matlab datenum)',t0) nc_attput(ROMS_frc_file,nc_global,'ini date (ROMS daynumber)',t0-datenum(1968,5,23)) nc_attput(ROMS_frc_file,nc_global,'components',ROMSnames_att) %dims dimname={'tide_period'; 'xi_rho'; 'eta_rho'; 'xi_u'; 'eta_u'; 'xi_v'; 'eta_v'}; dimlen=[Nharmonics; L+1; M+1; L; M+1; L+1; M]; for i=1:length(dimname); nc_adddim(ROMS_frc_file,dimname{i},dimlen(i)); end %vars varname={'tide_period'; 'tide_Ephase'; 'tide_Eamp';... 'tide_Cmin'; 'tide_Cmax'; 'tide_Cangle'; 'tide_Cphase';... 'tide_Uamp'; 'tide_Uphase'; 'tide_Vamp'; 'tide_Vphase'}; vartype={'double'; 'double'; 'double'; 'double'; 'double'; 'double';... 'double'; 'double'; 'double'; 'double'; 'double'}; vardims={{'tide_period'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'} {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}; {'tide_period','eta_rho','xi_rho'}}; nVar=length(varname); for i=1:nVar Var.Name=varname{i}; Var.Nctype=vartype{i}; Var.Dimension=vardims{i}; nc_addvar(ROMS_frc_file,Var); end %atts attnames={{'long_name','units'}; %tide_period {'long_name','units'}; %tide_Ephase {'long_name','units'}; %tide_Eamp {'long_name','units'}; %tide_Cmin {'long_name','units'}; %tide_Cmax {'long_name','units'}; %tide_Cangle {'long_name','units'}; %tide_Cphase {'long_name','units'}; %tide_Uamp {'long_name','units'}; %tide_Uphase {'long_name','units'}; %tide_Vamp {'long_name','units'}}; %tide_Vphase attvals={{'Tide angular period','hours'}; %tide_period {'Tide elevation phase angle','degrees'}; %tide_Ephase {'Tide elevation amplitude','meters'}; %tide_Eamp {'Tidal current ellipse semi-minor axis','meter second-1'}; %tide_Cmin {'Tidal current ellipse semi-major axis','meter second-1'}; %tide_Cmax {'Tidal current ellipse inclination angle','degrees between semi-major axis and east'}; %tide_Cangle {'Tidal current phase angle','degrees'}; %tide_Cphase {'Tidal current U-component amplitude','meters'} %tide_Uamp {'Tidal current U-component phase','degrees'} %tide_Uphase {'Tidal current V-component amplitude','meters'} %tide_Vamp {'Tidal current V-component phase','degrees'}}; %tide_Vphase [nAtt,k]=size(attnames); for n=1:nVar [k itot]=size(attnames{n}); for i=1:itot attname=char(attnames{n}); attname=strtrim(attname(i,:)); attval=attvals{n}; if ischar(attval{i}) attval=strtrim(char(attval{i})); %if character string else attval=attval{i}; %if numeric end nc_attput(ROMS_frc_file,varname{n},attname,attval) end end %Write to output file %"PRESERVE_FVD" suppresses the permute that nc_varput_tmw would do. %start/count/stride has been disabled in nc_varput_tmw. nc_varput(ROMS_frc_file,'tide_period',ROMSperiods,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Eamp',zamp,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Ephase',zpha,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Cmax',major,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Cmin',minor,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Cangle',inclination,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Cphase',phase,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Uamp',uamp,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Uphase',upha,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Vamp',vamp,'PRESERVE_FVD') nc_varput(ROMS_frc_file,'tide_Vphase',vpha,'PRESERVE_FVD') % Plot-check map of amplitude and ellipses figure Nplot=find(strcmp(MapPlotName,ROMSnames)); start=[Nplot-1 0 0]; count=[1 -1 -1]; zamp=nc_varget(ROMS_frc_file,'tide_Eamp',start,count); zpha=nc_varget(ROMS_frc_file,'tide_Ephase',start,count); Cmax=nc_varget(ROMS_frc_file,'tide_Cmax'); Cmax=squeeze(Cmax(Nplot,1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1)); Cmin=nc_varget(ROMS_frc_file,'tide_Cmin'); Cmin=squeeze(Cmin(Nplot,1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1)); Cang=nc_varget(ROMS_frc_file,'tide_Cangle'); Cang=squeeze(Cang(Nplot,1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1))*pi/180; cmpt=nc_attget(ROMS_frc_file,nc_global,'components'); %ROMS3: each comp. field has 4 places (3 in Roms_tools) zamp(maskR==0)=NaN; zpha(maskR==0)=NaN; contourf(lonR,latR,zamp) %must include rlon, rlat or ellipse won't work colorbar('EastOutside') slon=lonR(1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1); slat=latR(1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1); smask=maskR(1:ellipsePlotStride:M+1,1:ellipsePlotStride:L+1); tide_ellipse(ellipsePlotScale*smask.*Cmax,ellipsePlotScale*smask.*Cmin,smask.*Cang,slon,slat,'k'); title([char(ROMSnames(Nplot)) ' amplitude (m)']); lonmin=min(min(lonR)); lonmax=max(max(lonR)); latmin=min(min(latR)); latmax=max(max(latR)); axis([lonmin lonmax latmin latmax]) %Plot-check map of phase figure contourf(lonR,latR,zpha) title([char(ROMSnames(Nplot)) ' phase']); axis([lonmin lonmax latmin latmax]) colorbar('EastOutside') %Plot-check 28-day time series starting from t0. %Phase as required in the ROMS forcing file is equivalent to the sum %"g-(V+u)" where g is what some people call "greenwich phase lag". dist=spheric_dist(plot_timeseries_Lat,latR,plot_timeseries_Lon,lonR); [J,I]=find(dist==(min(min(dist)))); disp(['Predictions row,col: ' int2str(J) ' ' int2str(I)]); lon0=lonR(J,I); lat0=latR(J,I); disp(['Predictions coordinates ',num2str(lon0),' ',num2str(lat0)]) Tperiod=nc_varget(ROMS_frc_file,'tide_period'); zamp=nc_varget(ROMS_frc_file,'tide_Eamp'); zamp=squeeze(zamp(:,J,I)); zpha=nc_varget(ROMS_frc_file,'tide_Ephase'); zpha=squeeze(zpha(:,J,I)); t=datenum(t0:1/24:t0+28); %times (to plot) (GMT) fRPD=2*pi./(ROMSperiods/24); %frequency rad/day for i=1:Nharmonics phaselag=(pi/180)*zpha(i); phaselag=phaselag(:); ssh1(:,i)=zamp(i)*cos((t-t(1))*fRPD(i)-phaselag); end ssh=sum(ssh1,2); figure plot(t,ssh,'r') axis([t(1) t(end) -1.5 1.5]) datetick('x','ddmmm','keeplimits') ylabel('Sea level [m]') grid on hold on %t_xtide predictions if ~isempty(XtidePORT) xpred=t_xtide(XtidePORT,t+LTZ/24); %xtide constants are in local time zones (LTZ); plot in UT plot(t,xpred-mean(xpred),'-k') end %Observed time series in GMT if ~isempty(sealevel_obs_file) Observed=load(sealevel_obs_file); Obs.t=Observed(:,1); Obs.sl=Observed(:,2); rm=mean(Obs.sl(~isnan(Obs.sl))); it=find(Obs.t >=t0 & Obs.t<=t0+28); plot(Obs.t(it),Obs.sl(it)-rm,'-g') end