Fix signs + detailed comments

magic-sph · Jun 17, 2024 · a23c257 · a23c257
1 parent e9bb616
commit a23c257
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 4 deletions.
diff --git a/src/preCalculations.f90 b/src/preCalculations.f90
@@ -783,9 +783,10 @@ subroutine preCalc(tscheme)
          if ( amp_tide /= 0.0_cp ) then
             y20_norm = 0.5_cp  * sqrt(5.0_cp/pi)
             y22_norm = 0.25_cp * sqrt(7.5_cp/pi)
-            tide_fac20  = amp_tide / y20_norm * r_cmb**2
-            tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor
-            tide_fac22n = -7.0_cp * tide_fac22p                                ! Half factor carried over
+            tide_fac20  = amp_tide / y20_norm * r_cmb**2 ! (2,0,1) mode of Ogilvie 2014
+            tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor, (2,2,1) mode
+            tide_fac22n = -7.0_cp * tide_fac22p                                ! Half factor carried over, (2,2,3) mode,
+                                                                               ! has opposite sign to that of the other two (Polfliet & Smeyers, 1990)
          else
             tide_fac20  = 0.0_cp
             tide_fac22p = 0.0_cp

diff --git a/src/updateWP.f90 b/src/updateWP.f90
@@ -445,12 +445,17 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, &
                         end if
 
                         if ( amp_tide /= 0.0_cp ) then
+                           ! tide_fac22p -> (2,2,1)
+                           ! tide_fac22n -> (2,2,3)
+                           ! (2,2,3) must have the same signed frequency
+                           ! as (2,0,1) above while (2,2,1) has the opposite
+                           ! sign in the rotating frame (Ogilvie, 2014)
                            rhs1(1,2*lm-1,threadid)=rhs1(1,2*lm-1,threadid)        &
                            &                 + (tide_fac22p + tide_fac22n)        &
                            &                      /real(l1*(l1+1),kind=cp)        &
                            &                      * cos(omega_tide*(time))
                            rhs1(1,2*lm,threadid)  =rhs1(1,2*lm,threadid)          &
-                           &                 +  (tide_fac22p - tide_fac22n)       &
+                           &                 +  (-tide_fac22p + tide_fac22n)       &
                            &                       /real(l1*(l1+1),kind=cp)       &
                            &                       * sin(omega_tide*(time))
                         end if