Possible bug in set_tides.F

[lmp4]

I think I've found a bug in how the contribution to sea surface height is added from a tidal forcing file.

According to https://www.myroms.org/wiki/index.php/tide_start; the following formula is used to add the tides onto zeta.



Now looking at the code in set_tides.F;

Code: Select all

!-----------------------------------------------------------------------
!  Add tidal elevation (m) to sea surface height climatology.
!-----------------------------------------------------------------------
!
        Etide(:,:)=0.0_r8
        cff=2.0_r8*pi*(time(ng)-tide_start*day2sec)
        DO itide=1,NTC
          IF (Tperiod(itide).gt.0.0_r8) THEN
            omega=cff/Tperiod(itide)
            DO j=JstrR,JendR
              DO i=IstrR,IendR
                Etide(i,j)=Etide(i,j)+                                  &
     &                     ramp*SSH_Tamp(i,j,itide)*                    &
     &                     COS(omega-SSH_Tphase(i,j,itide))
#  ifdef MASKING
                Etide(i,j)=Etide(i,j)*rmask(i,j)
#  endif
              END DO
            END DO
          END IF
        END DO

Here there seems to be a discrepancy between the equation and the code as this is what I think the code is doing;

Code Euqation.png (2.27 KiB) Viewed 4127 times

So should the code be;

Code: Select all

   &                     COS(omega-(2*pi*SSH_Tphase(i,j,itide)/Tperiod(itide))

or is there another reason for leaving out the 2*pi/Ti factor ?

The reason why I'm so interested in the tides is because I want to add the tidal model into the observations so I can assimilate SSH whilst still including tides! Computing the tidal contribution to zeta from the tidal forcing harmonics in the same way as ROMS seemed like the best way to add the tides to the AVISO SSH observations for assimilating.
In doing so I wrote a matlab script to plot the tidal contribution to sea level that I'd like to partiality share;

Code: Select all

tide_start=1*day2sec; %tide_start needs to be in seconds among all the time variables including Tperiod

Tperiod=nc_read(TideFile,'tide_period');

SSH_Tphase=nc_read(TideFile,'tide_Ephase')/360; 

SSH_Tamp=nc_read(TideFile,'tide_Eamp'); Etide=zeros(iend,jend);
    
    cff=2*pi*(Time(n)-(tide_start));
    
    for itide=1:NTC
        if Tperiod(itide)>0
              omega=cff/(Tperiod(itide)*hour2sec);
            for j=1:jend
                for i=1:iend
&                             Etide(i,j)=Etide(i,j)+(SSH_Tamp(i,j,itide)*cos(omega-(SSH_Tphase(i,j,itide))));  % OLD VERSION!
                             Etide(i,j)=Etide(i,j)+(SSH_Tamp(i,j,itide)*cos(omega-((2*pi*SSH_Tphase(i,j,itide))/(Tperiod(itide)*hour2sec)))) % CORRECTED VERSION!;

                    Etide(i,j)=Etide(i,j)*mask(i,j);
                end
            end
        end
    end

The following figure, computed using the matlab script above, shows the tidal contribution to zeta in meters before (left) and after (right) the factor has been applied. Note the field has become smooth in the corrected version instead of disjointed as was the case for all different times I calculated.

I've also tried looking at other grids and other tidal forcing but the changes from 'corrected' to the current code dramatically vary from different locations Some were more simply just different and didn't experience this disjointed to smooth transition which makes me second guess my logic I am still quite the amateur at ocean modeling so if I have made a mistake that's quite obvious I won't feel at all bad to be corrected! Thanks!

[simion1232006]

The code should be correct. If my understanding is right, the phase angle should not have an a factor of period in it. Just make sure you have the correct phase angle input, and I think you will be fine.

[lmp4]

If the ROMS code is correct, then is the formula presented in https://www.myroms.org/wiki/index.php/tide_start wrong?
My phase angel is in degrees units from input so I assume somewhere within the code other from here it is converted into (degree/360) the standard units used as stated in mod_tides.F. In any case if the ROMS code is correct then my matlab code / tidal forcing must be wrong as I can't seem to avoid the large discrepancies in my zeta calculation.

[lmp4]

Ok I've made a vital mistake! I didn't think to double check the units in varinfo.dat instead I just looked at the information given in mod_tides.F....

Code: Select all

'tide_Ephase'                                      ! Input
  'tidal elevation phase angle'
  'degrees'                                        ! [radians]
  'SSH_Tphase, scalar, series'
  'tide_period'
  'idTzph'
  'r2dvar'
  0.017453292519943295                             ! pi/180

So when reverting to the original ROMS code and simply multiplying my SSH_Ehpase by (pi/180) such as in varinfo.dat, instead of (1/360) as I was doing, my output returns to a smooth field!

Tides_correct_units.png (8.77 KiB) Viewed 4041 times

However this only makes the choice between the two even harder as the 2*pi/Tperiod(itide) correction is also a smooth field just a different shape.... The 24hr variability in the original code version does seem more correct than my proposed change but I still cannot wrap my head around the difference between formula and the code =(