program TSUCAST c GOAL c This program is designed to generate forecasts of later tsunami waves -- c specifically, water-level curves encompassing the highest and lowest water c levels for the remainder of the tsunami. The forecasts are based on data c from real-time reporting tide gages and are for the regions immediately c adjacent to the respective gages. c The program is meant to be converted into a module within the larger data c processing/data display systems (DPDDS) at the NOAA and other tsunami c warning centers. These systems provide the tide gage observations and tidal c predictions to the TSUCAST module and display these quantities and the c forecasts in a convenient format. c [Reference: Mofjeld, H.O., F.I. Gonzalez, E.N. Bernard and J.C. Newman, c Forecasting the heights of later waves in Pacific-wide tsunamis, c Natural Hazards, in press] c FORTRAN 77 CODE c This is a stand-alone FORTRAN 77 program. It inputs a pair of time series c from a given tide gage: the water level data 'dat' & the corresponding c tidal predictions 'tide'. Test datasets and an archival version of this c program are available at PMEL's anonymous ftp site. Instructions for c downloading these files are given later in this documentation (section c TEST CASES). c TSUCAST uses some VAX/VMS statement conventions. c (Use the appropriate compiler options for non-VMS operating systems) c INPUT/OUTPUT c To obtain forecasts for different tide gages and tsunamis, the I/O c filenames are changed and rerunning the program. This form makes it c easy to test TSUCAST after it is integrated into the larger DPDDS. c The I/O is done in the subroutines 'input' and 'output'. The filenames c are specified externally. Within the program, these are given generic c filenames ('datafile', 'tidefile' and 'outfile') in the 'open' statements. c TIME BASE c It is assumed that absolute time stamps will be imposed on the forecast c values by the DPDDS. Internally, TSUCAST uses times (in hours) relative to c the time of the first data value. c Different tide gage systems and reporting schemes may have different c sampling intervals. Hence, the sampling interval 'dt' needs to be known c and given to the program. In TSUCAST, 'dt' (in minutes) is read with an c 'accept' statement; 'dt' is then converted internally to 'dthr' (in hours). c The input sampling interval is specified in minutes because this is easier c to understand and less prone to errors. The wave height forecasts are made c at the same sampling interval. c TIDAL PREDICTIONS c Predicted tides are needed by TSUCAST for two reasons. The first is to remove c the tides from the 2-4 hr data segment that is used to compute the tsunami- c band variance and, hence, the amplitudes of wave height envelopes. The c second reason is to include the tides in the forecasts. The National Ocean c Service displays predicted and observed tides for real-time reporting tide c gages (http://co-ops.nos.noaa.gov). The tidal prediction software can also c be resident at the warning centers. c TSUCAST needs separate tidal predictions for each tide gage. The predictions c must be for at least 52 hrs after the onset of the tsunami at the gage. c DETAILS OF THE METHOD c The forecasts begin 4 hrs after the tsunami onset 't0' and are for 48 hrs c (i.e., forecast ends at time=t0+52 hrs). c The onset 't0' is detected externally (e.g., by a tsunami trigger on the c tide gage or by the travel time of the tsunami to the gage). c Extreme wave envelopes (enclosing the highest & lowest waves) are added c to the predicted tide to get the total forecast. c The water level units are the same as the tide gage data and the predicted c tides (all the water level information must have consistent units and have c the same reference level; the forecasts are relative to the same datum). c The forecasts are made under the assumption that the tsunami energy decays c exponentially with time with an e-folding time tau of 48 hours. c forep(k)=tide(k)+A*sigma*exp(-t/tau) for highest wave peaks (high waters) c & foret(k)=tide(k)-A*sigma*exp(-t/tau) for deepest wave troughs (low waters) c where A=3.0 and t=0 at time=t0+4 hr (beginning of forecast period). c Here, the standard deviation 'sigma' is computed over the data interval c DEL=[t0+2,t0+4] c and is found from the difference between the data and tidal predictions c diff(k)=dat(k)-tide(k) for each k in the data interval DEL c sum=SUM[ diff(k)], where the sum is over the data in DEL c sigma=sqrt[sum/N], where N is the number of data values in DEL c Note that the first two hours of tsunami [t0,t+2) is not used in the c computation of 'sigma'. c It is assumed that the DPDDS scansthe water level data for data gaps c before launching TSUCAST. If a gap exists in the t=2-4 hr data window, c the process can be delayed until 2 hrs of continuous data is available c from the gage. c TEST CASES c Test data to check the code are available at the PMEL ftp site as ASCII c files. The sample output (data file and Postscript graph) is also available c from the same site, as is a copy of this program. One test case is based c on the 1952 Kamchatka tsunami measured (0.6 min sampling) by the tide gage c at Crescent City, California. The second case is based on theoretical c sinusoidal waves (1.0 min sampling). c The test cases can be run after downloading the input data files and c sample output file via anonymous ftp. This is done using the following c responses (case sensitive) to the ftp prompts: c ftp ftp.pmel.noaa.gov c login [use only if the first ftp prompt isn't for Name or username] c anonymous c [use your full e-mail address for the password] c cd tsunami/tsucast/ c mget * c quit c To compile, link and run using VAX VMS( 1952 Kamchatka tsunami): c for tsucast c link tsucast c assign cc_52_obs.dat datafile c assign cc_52_tides.dat tidefile c assign cc_52_forecast.dat outfile c run tsucast c 0.6 c [Outfile file is named 'cc_52_forecast.dat'] c To compile and run using SunOS (Unix): c f77 -o tsucast.e -xl tsucast.f c ln -s cc_52_obs.dat datafile c ln -s cc_52_tides.dat tidefile c tsucast.e c 0.6 c mv outfile cc_52_forecast.dat c [Output filename is then converted to 'cc_52_forecast.dat'] c For the sinusoidal case, substitute 'sinusoidal_test' for 'cc_52' c in the filenames and the time sampling parameter '1.0' for '0.6' . c CREDITS c The development of the forecast method was funded in part by the c U.S. National Tsunami Hazard Mitigation Program, the Deputy Under c Secretary of Defense for Space Integration (in support of the Pacific c Disaster Center via funds administered through the NASA NRA-98-OES-13 c Research Announcement), and the NOAA/Pacific Marine Environmental c Laboratory. c Code written by H.O. Mofjeld NOAA/PMEL/Tsunami Research Program c E-mail: mofjeld@pmel.noaa.gov c Version 1.0 (written June 12, 2000) c DECLARATIONS integer NMAX,KK2 parameter (NMAX=12500,KK2=1000) integer k,k1,k2,kmax,N real*4 dat(NMAX),tide(NMAX),diff(KK2) real*4 a,tau real*4 sum,sigma,envp,envt real*4 dt,t,expon(NMAX) real*4 forep(NMAX),foret(NMAX) data a /3.0/, tau /48.0/ c SAMPLING INTERVAL c Inputting the sampling interval (this flexibility is provided in case c different gages, reporting systems, and/or historical data sets have c different sampling intervals for the water level data) c The units of the sampling interval are minutes print *, ' Please input sampling interval (min, real) ' print *, ' (0.25 min minimum; 2.0 min maximum) ' accept *, dt if((dt.ge.0.25).and.(dt.le.2.0)) then dthr=dt/60.0 ! converting to time (hr) else print *, ' Sampling interval dt (min) is outside [0.25,2.0] ' print *, ' Program TSUCAST terminated ' stop ! Abnormal termination endif k1=1.5+2.0/dthr ! computing array limits based on sampling interval k2=1.5+4.0/dthr N=k2-k1 kmax=1.5+52.0/dthr if(kmax.gt.NMAX) then kmax=NMAX ! safeguarding length of time series print *, ' series length kmax too large, reset to NMAX ' else endif c EXPONENTIAL DECAY ARRAY c Computing the unit-amplitude exponential envelope for the prototype c (the envelope is same for every tide station having the same sampling c interval; operationally, this array could be pre-computed) do k=k2,kmax t=dthr*(k-k2) expon(k)=exp(-t/tau) end do ! k c DATA INPUT c (filenames supplied externally) call input(k2,kmax,dat,tide) c COMPUTING THE EXTREME WAVE FORECASTS do k=k1,k2-1 ! computing the difference series diff(k)=dat(k)-tide(k) end do ! k do k=k1,k2-1 ! computing the sum sum=sum+diff(k)**2 end do ! k if(sum.gt.0.0) then ! computing sigma sigma=sqrt(sum/(N+0.0)) else print *, ' Summed squares of differences ' print *, ' not positive definite ' print *, ' (check input data) ' print *, ' Program TSUCAST terminated ' stop endif envp= a*sigma ! computing the envelope coefficients envt=-a*sigma do k=k2,kmax ! computing the forecasts (arrays) forep(k)=tide(k)+envp*expon(k) foret(k)=tide(k)+envt*expon(k) end do ! k c OUTPUT call output(kmax,dat,tide,forep,foret) print *, ' Program TSUCAST completed ' stop ! Normal termination of program forecast end ! tsucast subroutine input(k2,kmax,dat,tide) c Inputting 'dat' and 'tide' arrays c Test case c This version of the input subroutine is set up for a format c that reads pre-existing tide gage data of historical tsunamis c that are available at PMEL c A separate program was required to pre-compute the tidal c predictions 'tide' for other cases. integer k,k2,kmax real*4 dat(1),tide(1) open(unit=5,file='datafile',readonly,status='old') read (5,*,err=10,end=20) ! skipping 3-line header read (5,*,err=10,end=20) read (5,*,err=10,end=20) read (5,1,err=10,end=20) (dat(k),k=1,k2-1) 1 format(f10.4) close(unit=5) open(unit=5,file='tidefile',readonly,status='old') read (5,*,err=30,end=40) ! skipping 3-line header read (5,*,err=30,end=40) read (5,*,err=30,end=40) read (5,1,err=30,end=40) (tide(k),k=1,kmax) close(unit=5) return ! Normal return c Problems encountered while inputting reading input files 10 print *, ' Error encountered while reading datafile ' go to 50 20 print *, ' End of file encountered while reading datafile ' go to 50 30 print *, ' Error encountered while reading tidefile ' go to 50 40 print *, ' End of file encountered while reading tidefile ' 50 print *, ' Program TSUCAST terminated ' stop ! Abnormal stop while inputting data or tide arrays end ! input subroutine output(kmax,dat,tide,forep,foret) c Outputting the forecasts (including the 4 hrs before the forecast c begins). Only the first 4 hrs of the observations are valid in the c output file, since these data are all that are available when the c forecast is made. However, the predicted tides are available for c the full time period (t=0-52 hrs). integer k,kmax real*4 dat(1),tide(1),forep(1),foret(1) open(unit=7,file='outfile',status='new') write (7,2) (k,dat(k),tide(k),forep(k),foret(k),k=1,kmax) 2 format( i10,4f10.4) close(unit=7) print *, ' Forecast outputted successfully ' return end ! output