tidal analysis of large output files

 Posts: 20
 Joined: Fri Dec 16, 2011 3:14 pm
 Location: School of Ocean Sciences
tidal analysis of large output files
Hello all,
I have been using ROMS for (barotropic) tidal simulations. For high resolution tidal models, big data files should be first transferred from the computer cluster (which does not have MATLAB due to licence issue) to PC, and then locally analysed. This task is usually time consuming and causes some problems, if one tries to do the tidal analysis for the whole domain.
I was wondering, if there is a more convenient way to do this task in ROMS (sorry, if this topic has been already addressed somewhere that I couldn't find). Postprocessing of the ROMS outputs using FORTRAN would be an option, which I would consider if there was not a better way. Also, I have seen in other hydrodynamic models that the original code has an option for the tidal analysis. I was wondering if some day this would be added to ROMS?
I appreciate your comments
Thanks
I have been using ROMS for (barotropic) tidal simulations. For high resolution tidal models, big data files should be first transferred from the computer cluster (which does not have MATLAB due to licence issue) to PC, and then locally analysed. This task is usually time consuming and causes some problems, if one tries to do the tidal analysis for the whole domain.
I was wondering, if there is a more convenient way to do this task in ROMS (sorry, if this topic has been already addressed somewhere that I couldn't find). Postprocessing of the ROMS outputs using FORTRAN would be an option, which I would consider if there was not a better way. Also, I have seen in other hydrodynamic models that the original code has an option for the tidal analysis. I was wondering if some day this would be added to ROMS?
I appreciate your comments
Thanks
Re: tidal analysis of large output files
Probably best to do the tidal analysis using FORTRAN code on the supercomputing cluster. I did this when I analysed an 8 Terabyte dataset for 3D global tides from HYCOM. Foreman's FORTRAN harmonic analysis code uses the same algorithm as matlab t_tide but doesn't provide 95% confidence intervals.

 Posts: 20
 Joined: Fri Dec 16, 2011 3:14 pm
 Location: School of Ocean Sciences
Re: tidal analysis of large output files
Thanks for the reply. As I mentioned, I was wondering about something more convenient than extracting the data from a netcdf file and then doing the tidal analysis both in FORTRAN.
Re: tidal analysis of large output files
It think it would be rather difficult to implement harmonic tidal analysis within the code since it would require that they entire time series be held in memory. For a large domain that could impose some serious memory usage issues; especially if one wanted to study the 3D tidal currents. You'd probably want to interpolate the 3D currents into zspace prior to the tidal analysis.
I think it's easiest to write a netcdf frontend for Foreman's classical harmonic tidal analysis code (or his versatile tidal analysis code) and postprocess the model output. I'm sure there are other FORTRAN codes out there that could be used as well. When I was doing my work on HYCOM I wrote a frontend for Foreman's code so that is could input FORTRAN binaries from the HYCOM output, removed a few legacy FORTRAN 77 statements, and implemented OMP. I ported everything to Texas and ran it there. It executed much faster than matlab and all I had to do was download the output from the harmonic analysis.
I think it's easiest to write a netcdf frontend for Foreman's classical harmonic tidal analysis code (or his versatile tidal analysis code) and postprocess the model output. I'm sure there are other FORTRAN codes out there that could be used as well. When I was doing my work on HYCOM I wrote a frontend for Foreman's code so that is could input FORTRAN binaries from the HYCOM output, removed a few legacy FORTRAN 77 statements, and implemented OMP. I ported everything to Texas and ran it there. It executed much faster than matlab and all I had to do was download the output from the harmonic analysis.
Re: tidal analysis of large output files
Hi Hashemi,
If you are comfortable modifying the ROMS code, it is not too difficult to do the harmonic analysis during the run. The steps are moreorless as follows:
 Select a set of tidal frequencies for analysis
 Determine the duration of run needed to separate the selected frequencies
 Start the model run, and wait for spin up to finish
 After spin up, accumulate sums of the desired fields times sines and cosines at the selected frequencies
 At the end of the run, solve a small linear algebra problem to convert the accumulated sums into inphase and quadrature harmonic constants
 Finally, some minor gymnastics might be needed to rotate the phases to the correct reference
The most significant computational demand is likely to be the additional memory to store the accumulations of the desired fields. Feel free to contact me off list, ezaron@pdx.edu, if you need help with the implementation.
Ed
If you are comfortable modifying the ROMS code, it is not too difficult to do the harmonic analysis during the run. The steps are moreorless as follows:
 Select a set of tidal frequencies for analysis
 Determine the duration of run needed to separate the selected frequencies
 Start the model run, and wait for spin up to finish
 After spin up, accumulate sums of the desired fields times sines and cosines at the selected frequencies
 At the end of the run, solve a small linear algebra problem to convert the accumulated sums into inphase and quadrature harmonic constants
 Finally, some minor gymnastics might be needed to rotate the phases to the correct reference
The most significant computational demand is likely to be the additional memory to store the accumulations of the desired fields. Feel free to contact me off list, ezaron@pdx.edu, if you need help with the implementation.
Ed
Re: tidal analysis of large output files
I think a better solution for the issue is using either python or ncl or maybe perl, which you can use it in cluster without violating licence. The desired code is available in pyroms and hongkong workshop.
Re: tidal analysis of large output files
Thanks for the additional input. My biggest concern with implementing tidal analysis in ROMS is the memory requirement. If you wanted to include the four largest semidiurnal and four largest diurnal constituents (M2, S2, N2, K2, K1, O1, P1, and Q1) you would need 183 days of hourly data. That would require a large amount of memory for a large domain. Would work for selected locations or for shorter runs with fewer tidal constituents, however. Using python, ncl, or perl is certainly a good option. Be easier than writing a netcdf frontend in FORTRAN in order to read the ROMS output files.
Re: tidal analysis of large output files
The AVERAGES_DETIDE option does tidal harmonic analysis on the fly while the model runs.
It is not necessary to hold the entire time series in memory at all.
Harmonic analysis needs only the sum of the product of each time series with sin(wt) and cos(wt) for every harmonic. AVERAGES_DETIDE accumulates those summations as the model runs and on each averaging interval the Normal Equations are solved to complete the harmonic fit. With AVERAGES_DETIDE the internal harmonic fit is used to output detided averages which can avoid some problems with tidal aliasing in simple boxcar averages.
See this post in the source code trac:
https://www.myroms.org/projects/src/ticket/115
The new detided output is written to variables zeta_detided, ubar_detided etc. in the averages file as controlled by the flags in ocean.in:
Aout(idu3dD) == T ! detided 3D Uvelocity
Aout(idv3dD) == T ! detided 3D Vvelocity
Aout(idu2dD) == T ! detided 2D Uvelocity
Aout(idv2dD) == T ! detided 2D Vvelocity
Aout(idFsuD) == T ! detided freesurface
The harmonic fit coefficients themselves are appended to the tide forcing file as variables zeta_tide, ubar_tide etc. including the 3D variables u_tide, temp_tide etc. There you will also see terms like:
double CosWCosW(tide_period, tide_period) ;
CosWCosW:long_name = "timeaccumulated COS(omega(k)*t)*COS(omega(l)*t) matrix" ;
CosWCosW:units = "radians" ;
which are the accumulated products. The products and the harmonics are continually updated as the model runs.
The longer the model runs the more stable and accurate become the fitted harmonics.
There is a substantial CPU demand with running the AVERAGES_DETIDE option, so if you use this option you'll probably only want to run as long as you need to to separate out the harmonics you care about. Not selecting the 3D variables for output will lessen the computation.
It is not necessary to hold the entire time series in memory at all.
Harmonic analysis needs only the sum of the product of each time series with sin(wt) and cos(wt) for every harmonic. AVERAGES_DETIDE accumulates those summations as the model runs and on each averaging interval the Normal Equations are solved to complete the harmonic fit. With AVERAGES_DETIDE the internal harmonic fit is used to output detided averages which can avoid some problems with tidal aliasing in simple boxcar averages.
See this post in the source code trac:
https://www.myroms.org/projects/src/ticket/115
The new detided output is written to variables zeta_detided, ubar_detided etc. in the averages file as controlled by the flags in ocean.in:
Aout(idu3dD) == T ! detided 3D Uvelocity
Aout(idv3dD) == T ! detided 3D Vvelocity
Aout(idu2dD) == T ! detided 2D Uvelocity
Aout(idv2dD) == T ! detided 2D Vvelocity
Aout(idFsuD) == T ! detided freesurface
The harmonic fit coefficients themselves are appended to the tide forcing file as variables zeta_tide, ubar_tide etc. including the 3D variables u_tide, temp_tide etc. There you will also see terms like:
double CosWCosW(tide_period, tide_period) ;
CosWCosW:long_name = "timeaccumulated COS(omega(k)*t)*COS(omega(l)*t) matrix" ;
CosWCosW:units = "radians" ;
which are the accumulated products. The products and the harmonics are continually updated as the model runs.
The longer the model runs the more stable and accurate become the fitted harmonics.
There is a substantial CPU demand with running the AVERAGES_DETIDE option, so if you use this option you'll probably only want to run as long as you need to to separate out the harmonics you care about. Not selecting the 3D variables for output will lessen the computation.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
Re: tidal analysis of large output files
Thanks for correcting me. Now that I've thought about it more you don't need to hold the time series in memory. I'll take a look at the AVERAGES_DETIDE option.

 Posts: 20
 Joined: Fri Dec 16, 2011 3:14 pm
 Location: School of Ocean Sciences
Re: tidal analysis of large output files
Thank you all for your replies and help. AVERAGES_DETIDE or using python seems to solve the problem. With regard to reading all data into memory, I think even in the MATLAB tidal analysis we can avoid that. I mean, we can extract the time series at a single node and store it in a dummy variable and after computing harmonic constants move to another node. Although, based my experience, this approach is slower, it doesn't need a huge memory to store the whole hydrodynamic field in the memory.
 jivica
 Posts: 130
 Joined: Mon May 05, 2003 2:41 pm
 Location: The University of Western Australia, Perth, Australia
Re: tidal analysis of large output files
If you decide to go the way of python or matlab then you can use pp (parallel python) or parallel toolbox in matlab and distribute the node_loop along the cluster (cores); doing that you are minimizing memory footprint and speedup calcs. Note that MPI is reducing memory footprint and not openmp (cores at the one mbo). Actually, when I think it over, it is easy to write simple parallel fortran prog (or adjust existing one) to read ROMS netCDF and use MPI_SEND, _RECEIVE, there are numerous examples in the distribution of openMPI...
Cheers,
Ivica
Cheers,
Ivica

 Posts: 82
 Joined: Mon Aug 16, 2004 8:47 pm
 Location: U.S. Geological Survey, Woods Hole
 Contact:
Re: tidal analysis of large output files
I tried AVERAGES_DETIDE, and the averages file contained "zeta_detided" etc but none of the mentioned harmonic fit coefficients. Did I mess something up with the file writing specs (relevant ones below)? I can provide more info but assumed this is a straightforward error on my part.
DT = 5.0
NTIMES == 362880
NTSAVG == 120900
NAVG == 241920
NDEFAVG == 0
NTSDIA == 1
NDIA == 720
NDEFDIA == 0
DT = 5.0
NTIMES == 362880
NTSAVG == 120900
NAVG == 241920
NDEFAVG == 0
NTSDIA == 1
NDIA == 720
NDEFDIA == 0

 Posts: 82
 Joined: Mon Aug 16, 2004 8:47 pm
 Location: U.S. Geological Survey, Woods Hole
 Contact:
Re: tidal analysis of large output files
ignore that...just reread Wilkin's post more carefully.
Re: tidal analysis of large output files
Looking at the output from a simulation run with the AVERAGES_DETIDE option, I see:
so there are 9 tide_periods, and 19 harmonics. I'm guessing that these 19 harmonics are the 9 tide_periods + 9 sum frequencies (0.5 tide_periods) + mean. I say I'm guessing because there are no variables called "harmonic_period" or "harmonic_names" or similar. The only variables that depend on the dimension "harmonics" are the analyzed variables zeta_tide, ubar_tide, etc.
How do I determine what the frequencies of the 19 harmonics are?
Thanks,
Rich
Code: Select all
$ncdump h bbleh_base_detide.nc
netcdf bbleh_base_detide {
dimensions:
xi_rho = 160 ;
eta_rho = 800 ;
tide_period = 9 ;
eta_u = 800 ;
eta_v = 799 ;
xi_v = 160 ;
xi_u = 159 ;
s_rho = 7 ;
harmonics = 19 ;
variables:
double tide_period(tide_period) ;
tide_period:long = "tide angular period" ;
tide_period:units = "hours" ;
tide_period:field = "tide_period, scalar" ;
...
double zeta_tide(harmonics, eta_rho, xi_rho) ;
zeta_tide:long_name = "timeaccumulated freesurface tide harmonics" ;
zeta_tide:units = "meter" ;
How do I determine what the frequencies of the 19 harmonics are?
Thanks,
Rich
Re: tidal analysis of large output files
The comments in the preamble in Modules/mod_tides.F explains the order of the harmonics data is mean followed by sin(wt) terms then cos(wt) terms. Within each group the order of the omega (w) harmonics is as given by the user in the tides forcing file.
Code: Select all
# if defined AVERAGES_DETIDE && defined AVERAGES
! !
! Detided timeaveraged fields via leastsquares fitting. Notice that !
! the harmonics for the state variable have an extra dimension of !
! size (0:2*NTC) to store several terms: !
! !
! index 0 mean term (accumulated sum) !
! 1:NTC accumulated sine terms !
! NTC+1:2*NTC accumulated cosine terms !
! !
! CosW_avg Current timeaverage window COS(omega(k)*t). !
! CosW_sum Timeaccumulated COS(omega(k)*t). !
! SinW_avg Current timeaverage window SIN(omega(k)*t). !
! SinW_sum Timeaccumulated SIN(omega(k)*t). !
! CosWCosW Timeaccumulated COS(omega(k)*t)*COS(omega(l)*t). !
! SinWSinW Timeaccumulated SIN(omega(k)*t)*SIN(omega(l)*t). !
! SinWCosW Timeaccumulated SIN(omega(k)*t)*COS(omega(l)*t). !
! !
! ubar_detided Timeaveraged and detided 2D umomentum. !
! ubar_tide Timeaccumulated 2D umomentum tide harmonics. !
! vbar_detided Timeaveraged and detided 2D vmomentum. !
! vbar_tide Timeaccumulated 2D vmomentum tide harmonics. !
! zeta_detided Timeaveraged and detided freesurface. !
! zeta_tide Timeaccumulated freesurface tide harmonics. !
# ifdef SOLVE3D
! t_detided Timeaveraged and detided tracers (T,S). !
! t_tide Timeaccumulated 3D tracers (T,S) tide harmonics. !
! u_detided Timeaveraged and detided 3D umomentum. !
! u_tide Timeaccumulated 3D umomentum tide harmonics. !
! v_detided Timeaveraged and detided 3D vmomentum. !
! v_tide Timeaccumulated 3D vmomentum tide harmonics. !
# endif
! !
# endif
!=======================================================================
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
Re: tidal analysis of large output files
Excellent. I forgot one of the fundamentals of modeling: when in doubt, look at the code!
So if we are forcing with M2 only, but want to look at nonlinear tidal effects, we need to add M4, M6, M8 to the forcing file so that they get included in the leastsquares fit, right?
So if we are forcing with M2 only, but want to look at nonlinear tidal effects, we need to add M4, M6, M8 to the forcing file so that they get included in the leastsquares fit, right?
Re: tidal analysis of large output files
I think you can just put these tidal constituents (M4,M6 and M8) into the tide forcing file but set their amplitude to zero.
Attached is a matlab code to read the output tidal forcing files and convert the harmonic summation terms to phase and amplitude for each variable.
And a time series of zeta is posted here. The tidal analysis runs for about 150 days with 8 constituents in the Chesapeake Bay.
Attached is a matlab code to read the output tidal forcing files and convert the harmonic summation terms to phase and amplitude for each variable.
And a time series of zeta is posted here. The tidal analysis runs for about 150 days with 8 constituents in the Chesapeake Bay.
 Attachments

 harm_roms.m
 (6.49 KiB) Downloaded 141 times
Re: tidal analysis of large output files
ROMS folk,
We tried using bzhang's harm_roms.m script to calculate the tidal amplitudes and phases from the AVERAGES_DETIDE output, but our amplitudes are more than a factor or 30 too high. bzhang couldn't see anything wrong. In case someone is willing to take a look, the output file is at:
http://geoport.whoi.edu/thredds/fileSer ... _detide.nc and our version of harm_roms.m to read it is attached below.
We ran the model from a cold start, and used NTSAVG to start writing to the AVERAGES file after an initial spinup period, figuring we didn't want that in the tidal analysis. We then saved the rest of the run into a single AVERAGES file.
Does this cause a problem for analyzing the AVERAGES_DETIDE output?
Here's specifically what we used:
We tried using bzhang's harm_roms.m script to calculate the tidal amplitudes and phases from the AVERAGES_DETIDE output, but our amplitudes are more than a factor or 30 too high. bzhang couldn't see anything wrong. In case someone is willing to take a look, the output file is at:
http://geoport.whoi.edu/thredds/fileSer ... _detide.nc and our version of harm_roms.m to read it is attached below.
We ran the model from a cold start, and used NTSAVG to start writing to the AVERAGES file after an initial spinup period, figuring we didn't want that in the tidal analysis. We then saved the rest of the run into a single AVERAGES file.
Does this cause a problem for analyzing the AVERAGES_DETIDE output?
Here's specifically what we used:
Code: Select all
! Input/Output parameters.
NRREC == 0
LcycleRST == T
NRST == 4500000
NSTA == 1
NFLT == 900
NINFO == 100
! Output history, average, diagnostic files parameters.
LDEFOUT == T
NHIS == 18000000
NDEFHIS == 1800
NTSAVG == 120900
NAVG == 241920
NDEFAVG == 0
NTSDIA == 1
NDIA == 720
NDEFDIA == 0
 Attachments

 harm_roms.m
 (6.16 KiB) Downloaded 137 times
 arango
 Site Admin
 Posts: 1102
 Joined: Wed Feb 26, 2003 4:41 pm
 Location: IMCS, Rutgers University
 Contact:
Re: tidal analysis of large output files
It is possible that your solution is not long enough to separate all the 9 tidal frequencies in your application. Recall that we have in ROMS the nonlinear interaction between all the frequencies and resonance. If you have an application with just M_{2} and S_{2}, it will take around 30 days of solution to separate those frequencies in your application. I expect that your application with 9 harmonics, it will take much much longer. Usually it will take months or couple of years of solution. ROMS was designed to improve the detide separation as the simulation becomes longer and longer. That's the reason why we have a restart capability for this.
In your case, the harmonics dimension (2*NCT+1) in the NetCDF file is 19 because you have NTC=9. The sequence is as follows:
The order of the harmonics is the same as specified in the input tidal forcing NetCDF file. In set_avg.F, we solve a leastsquares problem:
My tidal forcing NetCDF files always have the information about the frequencies in the global attributes:
so you know what is the order of frequencies inside ROMS.
In your case, the harmonics dimension (2*NCT+1) in the NetCDF file is 19 because you have NTC=9. The sequence is as follows:
Code: Select all
Harmonics(1) = mean
Harmonics(2:10) = sin(omega(k)*t)
Harmonics(11:19) = cos(omega(k)*t)
Code: Select all
! Compute detide leastsquares coefficients. Build coefficient squared
! matrix C(0:2*NTC,0:2*NTC) to invert. It is 2*NTC because we are
! solving for real components Ak and Bk. The zero rows and column is
! for the coefficients associated with the time mean.
!
! F(t) = Fmean + SUM [ Ak sin(omega(k)*t) ]
! + SUM [ Bk cos(omega(k)*t) ] for k=1:NTC
!
! In the code below, all the arrays are collapsed into a single
! dimension index such that:
!
! k=0 mean term
! k=1:NTC sine terms
! k=NTC+1:2*NTC cosine terms
Code: Select all
netcdf gom_tides_g {
dimensions:
xi_rho = 160 ;
eta_rho = 100 ;
tide_period = UNLIMITED ; // (4 currently)
...
// global attributes:
:type = "ROMS Forcing file" ;
:title = "ROMS Tidal Forcing from DATA/Model_tpxo7" ;
:grid_file = "/Users/arango/ocean/repository/Projects/gom/Data/g/gom_grid_g.nc" ;
:base_date = "days since 20000101 00:00:00" ;
:components = "S2, M2, K1, O1" ;
:source = "OTPS" ;
:history = "01Dec2013 20:06:24: Created by arango with write_tides.m" ;
data:
...
tide_period = 12.0000000048, 12.4206011981605, 23.934469605045,
25.819341694366 ;
Re: tidal analysis of large output files
I think the short simulation time may be the problem, as Hernan said. Attached are two figures with 60 days and 150 days of model runs with 8 tidal constituents in the CBOFS. The results are clearly different and 60 days are not enough to separate these tidal constituents. And the 60 days results look like that Rich described.
Re: tidal analysis of large output files
I follow wilkin's instruction to simulate tide only with M2 and S2 constitutes and the tide start time is set to 2000010100 ,then after three months running I use the harm_roms.m code provided by bzhang to analysis the expanded tide forcing file and got the tide_amp and tide_phase. but I don't know why the tide_phase values are all negative??
 Attachments

 cotide.jpg (29.92 KiB) Viewed 5560 times
Re: tidal analysis of large output files
Hi, WRH,
I see your simulated tidal phase is completely normal, the phase is out of atan2 in matlab and scaled by 180/pi. You just need to add a constant 360 degree to get a positive number. You can look at this paper for comparison: doi:10.1029/2002JC001451
I see your simulated tidal phase is completely normal, the phase is out of atan2 in matlab and scaled by 180/pi. You just need to add a constant 360 degree to get a positive number. You can look at this paper for comparison: doi:10.1029/2002JC001451
Re: tidal analysis of large output files
thank you for your reply!bzhang wrote:Hi, WRH,
I see your simulated tidal phase is completely normal, the phase is out of atan2 in matlab and scaled by 180/pi. You just need to add a constant 360 degree to get a positive number. You can look at this paper for comparison: doi:10.1029/2002JC001451