tidal analysis of large output files

[smrhashemi]

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

[ptimko]

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 3-D 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.

[smrhashemi]

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.

[ptimko]

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 3-D tidal currents. You'd probably want to interpolate the 3-D currents into z-space prior to the tidal analysis.

I think it's easiest to write a netcdf front-end for Foreman's classical harmonic tidal analysis code (or his versatile tidal analysis code) and post-process 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 front-end 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.

[ezaron]

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 more-or-less 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 in-phase 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

[p.heidary]

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 hong-kong workshop.

[ptimko]

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 semi-diurnal 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 front-end in FORTRAN in order to read the ROMS output files.

[wilkin]

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 U-velocity
Aout(idv3dD) == T ! detided 3D V-velocity
Aout(idu2dD) == T ! detided 2D U-velocity
Aout(idv2dD) == T ! detided 2D V-velocity
Aout(idFsuD) == T ! detided free-surface

The harmonic fit coefficients themselves are appended to the tide forcing file as variables zeta_tide, ubar_tide etc. including the 3-D variables u_tide, temp_tide etc. There you will also see terms like:

double CosWCosW(tide_period, tide_period) ;
CosWCosW:long_name = "time-accumulated 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 3-D variables for output will lessen the computation.

[ptimko]

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.

[smrhashemi]

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]

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

[nganju]

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

[nganju]

ignore that...just reread Wilkin's post more carefully.

[rsignell]

Looking at the output from a simulation run with the AVERAGES_DETIDE option, I see:

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 = "time-accumulated free-surface tide harmonics" ;
                zeta_tide:units = "meter" ;

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

[wilkin]

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 time-averaged fields via least-squares 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 time-average window COS(omega(k)*t).           !
!  CosW_sum     Time-accumulated COS(omega(k)*t).                      !
!  SinW_avg     Current time-average window SIN(omega(k)*t).           !
!  SinW_sum     Time-accumulated SIN(omega(k)*t).                      !
!  CosWCosW     Time-accumulated COS(omega(k)*t)*COS(omega(l)*t).      !
!  SinWSinW     Time-accumulated SIN(omega(k)*t)*SIN(omega(l)*t).      !
!  SinWCosW     Time-accumulated SIN(omega(k)*t)*COS(omega(l)*t).      !
!                                                                      !
!  ubar_detided Time-averaged and detided 2D u-momentum.               !
!  ubar_tide    Time-accumulated 2D u-momentum tide harmonics.         !
!  vbar_detided Time-averaged and detided 2D v-momentum.               !
!  vbar_tide    Time-accumulated 2D v-momentum tide harmonics.         !
!  zeta_detided Time-averaged and detided free-surface.                !
!  zeta_tide    Time-accumulated free-surface tide harmonics.          !
#  ifdef SOLVE3D
!  t_detided    Time-averaged and detided tracers (T,S).               !
!  t_tide       Time-accumulated 3D tracers (T,S) tide harmonics.      !
!  u_detided    Time-averaged and detided 3D u-momentum.               !
!  u_tide       Time-accumulated 3D u-momentum tide harmonics.         !
!  v_detided    Time-averaged and detided 3D v-momentum.               !
!  v_tide       Time-accumulated 3D v-momentum tide harmonics.         !
#  endif
!                                                                      !
# endif
!=======================================================================

[rsignell]

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 least-squares fit, right?

[bzhang]

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.

[rsignell]

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 spin-up 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

[arango]

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₂ and S₂, 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:

Code: Select all

Harmonics(1)     =  mean
Harmonics(2:10)  =  sin(omega(k)*t)
Harmonics(11:19) =  cos(omega(k)*t)

The order of the harmonics is the same as specified in the input tidal forcing NetCDF file. In set_avg.F, we solve a least-squares problem:

Code: Select all

!  Compute detide least-squares 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

My tidal forcing NetCDF files always have the information about the frequencies in the global attributes:

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 2000-01-01 00:00:00" ;
                :components = "S2, M2, K1, O1" ;
                :source = "OTPS" ;
                :history = "01-Dec-2013 20:06:24:  Created by arango with write_tides.m" ;

data:
...

 tide_period = 12.0000000048, 12.4206011981605, 23.934469605045, 
    25.819341694366 ;

so you know what is the order of frequencies inside ROMS.

[bzhang]

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.

[wrh]

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??

[bzhang]

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

[wrh]

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

thank you for your reply!