There are three steps, thanks to FERRET. (there would have been only one
step if I have to write all the code in C. But FERRET is great as you know.)

1. project winds into wave modes.
Under FERRET, run wave_comTaux_all.jnl. This program is only a 'main' program
which run a program called wave_comTaux_1mode.jnl. The latter does real dirty
calculation for each vertical mode. If you need only 4 vertical modes, run the
latter four times, as specified the first four lines of code.

2. integrate wind projection to oceanic response for each wave mode
Run simple cshell script run_all.csh, which execute kelvin.x and rossby.x.
Each execution takes an inputfile, which contains actual input data. In 
these inputfiles, you must specify the vertical wave number you want to
have manually. You may be able to make it more flexible.

In rossby.c and kelvin.c, you need specify totalSteps (time steps in your run),
damping coefficient (drag0, sorry for bad terminology). These two are in macro
list. You may change other macro values too. Also, there is no reason for 
integrating K-waves and R-waves separately. You can merge them if you like.

3. Compile model output into netCDF files
Run wave_long.jnl. This execution will generate four netCDF files each for
whole time series and seasonal cycle. The file names for the whole time series
are:
kelvin.cdf: forced K-waves;
kelvin_refr.cdf: reflected R-waves from eastern BD;
rossby.cdf: forced R-waves;
rossby_refk.cdf: reflected K-waves from western BD.

For seasonal cycle, file names are marked with 'clim'. In addition, you will
find R-waves files marked by '3mod' for '6mod', indicating the number of
meridional modes.

Finally, you can use these netCDF files. Since the reflection is computed assuming
infinite S-N straight boundary, you need use an efficiency parameter for actual
boundary reflection.

NOTE: you always need adjust time axis for your dataset. In FERRET code, adjust
def axis/t and set region. In C code, adjust the variable totalSteps.