Index: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90
===================================================================
--- NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90	(revision 12043)
+++ NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90	(revision 12045)
@@ -478,5 +478,5 @@
       sh_p = MOD( sh_p,  2*rpi )
 
-      cosI = 0.913694997 -0.035692561 *cos(sh_N)
+      cosI = 0.913694997_wp - 0.035692561_wp * COS( sh_N )
 
       sh_I = ACOS( cosI )
@@ -485,36 +485,33 @@
       sh_tgn2 = tan(sh_N/2.0)
 
-      at1=atan(1.01883*sh_tgn2)
-      at2=atan(0.64412*sh_tgn2)
-
-      sh_xi=-at1-at2+sh_N
-
-      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi
+      at1 = ATAN( 1.01883_wp * sh_tgn2 )
+      at2 = ATAN( 0.64412_wp * sh_tgn2 )
+
+      sh_xi = sh_N - at1 - at2
+
+      IF( sh_N > rpi ) sh_xi = sh_xi - 2.0_wp * rpi
 
       sh_nu = at1 - at2
 
-      !----------------------------------------------------------------------
-      ! For constituents l2 k1 k2
-      !----------------------------------------------------------------------
-
-      tgI2 = tan(sh_I/2.0)
-      P1   = sh_p-sh_xi
-
-      t2 = tgI2*tgI2
-      t4 = t2*t2
-      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
-
-      p = sin(2.0*P1)
-      q = 1.0/(6.0*t2)-cos(2.0*P1)
-      sh_R = atan(p/q)
-
-      p = sin(2.0*sh_I)*sin(sh_nu)
-      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
-      sh_nuprim = atan(p/q)
-
-      s2 = sin(sh_I)*sin(sh_I)
-      p  = s2*sin(2.0*sh_nu)
-      q  = s2*cos(2.0*sh_nu)+0.0727
-      sh_nusec = 0.5*atan(p/q)
+      ! For computation of tidal constituents L2 K1 K2
+      tgI2 = tan( sh_I / 2.0_wp )
+      P1   = sh_p - sh_xi
+      !
+      t2   = tgI2 * tgI2
+      t4   = t2 * t2
+      sh_x1ra = SQRT( 1.0 - 12.0 * t2 * COS( 2.0 * P1 ) + 36.0_wp * t4 )
+      !
+      p = SIN( 2.0_wp * P1 )
+      q = 1.0_wp / ( 6.0_wp * t2 ) - COS( 2.0_wp * P1 )
+      sh_R = ATAN( p / q )
+      !
+      p = SIN( 2.0_wp * sh_I ) * SIN( sh_nu )
+      q = SIN( 2.0_wp * sh_I ) * COS( sh_nu ) + 0.3347_wp
+      sh_nuprim = ATAN( p / q )
+      !
+      s2 = SIN( sh_I ) * SIN( sh_I )
+      p  = s2 * SIN( 2.0_wp * sh_nu )
+      q  = s2 * COS( 2.0_wp * sh_nu ) + 0.0727_wp
+      sh_nusec = 0.5_wp * ATAN( p / q )
       !
    END SUBROUTINE astronomic_angle
@@ -595,68 +592,79 @@
    RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
       !!----------------------------------------------------------------------
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kformula
-      !
-      REAL(wp)            :: zf
-      REAL(wp)            :: zs, zf1, zf2
-      CHARACTER(LEN=3)    :: clformula
+      !!                  ***  FUNCTION nodal_factort  ***
+      !!
+      !! ** Purpose : Compute amplitude correction factors due to nodal motion
+      !!
+      !! ** Reference :
+      !!       S58) Schureman, P. (1958): Manual of Harmonic Analysis and
+      !!            Prediction of Tides (Revised (1940) Edition (Reprinted 1958
+      !!            with corrections). Reprinted June 2001). U.S. Department of
+      !!            Commerce, Coast and Geodetic Survey Special Publication
+      !!            No. 98. Washington DC, United States Government Printing
+      !!            Office. 317 pp. DOI: 10.25607/OBP-155.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kformula
+      !
+      REAL(wp)            ::   zf
+      REAL(wp)            ::   zs, zf1, zf2
+      CHARACTER(LEN=3)    ::   clformula
       !!----------------------------------------------------------------------
       !
       SELECT CASE( kformula )
       !
-      CASE( 0 )                  !==  formule 0, solar waves
+      CASE( 0 )                  ! Formula 0, solar waves
          zf = 1.0
          !
-      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
-         zf=nodal_factort(78)
+      CASE( 1 )                  ! Formula 1, compound waves (78 x 78)
+         zf=nodal_factort( 78 )
          zf = zf * zf
          !
-      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235) 
-         zf1 = nodal_factort( 78)
+      CASE ( 4 )                 ! Formula 4,  compound waves (78 x 235) 
+         zf1 = nodal_factort( 78 )
          zf  = nodal_factort(235)
          zf  = zf1 * zf
          !
-      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
-         zf1 = nodal_factort(78)
+      CASE( 18 )                 ! Formula 18,  compound waves (78 x 78 x 78 )
+         zf1 = nodal_factort( 78 )
          zf  = zf1 * zf1 * zf1
          !
-      CASE( 20 )                 !==  formula 20, compound waves ( 78 x 78 x 78 x 78 )
-         zf1 = nodal_factort(78)
+      CASE( 20 )                 ! Formula 20, compound waves ( 78 x 78 x 78 x 78 )
+         zf1 = nodal_factort( 78 )
          zf  = zf1 * zf1 * zf1 * zf1
          !
-      CASE( 73 )                 !==  formule 73
-         zs = sin(sh_I)
-         zf = (2./3.-zs*zs)/0.5021
-         !
-      CASE( 74 )                 !==  formule 74
-         zs = sin(sh_I)
-         zf = zs * zs / 0.1578
-         !
-      CASE( 75 )                 !==  formule 75
-         zs = cos(sh_I/2)
-         zf = sin(sh_I) * zs * zs / 0.3800
-         !
-      CASE( 76 )                 !==  formule 76
-         zf = sin(2*sh_I) / 0.7214
-         !
-      CASE( 78 )                 !==  formule 78
-         zs = cos(sh_I/2)
-         zf = zs * zs * zs * zs / 0.9154
-         !
-      CASE( 149 )                !==  formule 149
-         zs = cos(sh_I/2)
-         zf = zs*zs*zs*zs*zs*zs / 0.8758
-         !
-      CASE( 215 )                !==  formule 215
-         zs = cos(sh_I/2)
-         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
-         !
-      CASE( 227 )                !==  formule 227 
-         zs = sin(2*sh_I)
-         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
-         !
-      CASE ( 235 )               !==  formule 235 
-         zs = sin(sh_I)
-         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
+      CASE( 73 )                 ! Formula 73 of S58
+         zs = SIN( sh_I )
+         zf = ( 2.0_wp / 3.0_wp - zs * zs ) / 0.5021_wp
+         !
+      CASE( 74 )                 ! Formula 74 of S58
+         zs = SIN(sh_I)
+         zf = zs * zs / 0.1578_wp
+         !
+      CASE( 75 )                 ! Formula 75 of S58
+         zs = COS( sh_I / 2.0_wp )
+         zf = SIN( sh_I ) * zs * zs / 0.3800_wp
+         !
+      CASE( 76 )                 ! Formula 76 of S58
+         zf = SIN( 2.0_wp * sh_I ) / 0.7214_wp
+         !
+      CASE( 78 )                 ! Formula 78 of S58
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs / 0.9154_wp
+         !
+      CASE( 149 )                ! Formula 149 of S58
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs * zs * zs / 0.8758_wp
+         !
+      CASE( 215 )                ! Formula 215 of S58 with typo correction (0.9154 instead of 0.9145)
+         zs = COS( sh_I/2 )
+         zf = zs * zs * zs * zs / 0.9154_wp * sh_x1ra
+         !
+      CASE( 227 )                ! Formula 227 of S58
+         zs = SIN( 2.0_wp * sh_I )
+         zf = SQRT( 0.8965_wp * zs * zs + 0.6001_wp * zs * COS( sh_nu ) + 0.1006_wp )
+         !
+      CASE ( 235 )               ! Formula 235 of S58
+         zs = SIN( sh_I )
+         zf = SQRT( 19.0444_wp * zs * zs * zs * zs + 2.7702_wp * zs * zs * cos( 2.0_wp * sh_nu ) + 0.0981_wp )
          !
       CASE DEFAULT
@@ -670,24 +678,25 @@
    FUNCTION dayjul( kyr, kmonth, kday )
       !!----------------------------------------------------------------------
-      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
-      !!----------------------------------------------------------------------
-      INTEGER,INTENT(in) ::   kyr, kmonth, kday
-      !
-      INTEGER,DIMENSION(12) ::  idayt, idays
-      INTEGER  ::   inc, ji
-      REAL(wp) ::   dayjul, zyq
-      !
-      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
-      !!----------------------------------------------------------------------
-      !
-      idays(1) = 0.
-      idays(2) = 31.
-      inc = 0.
-      zyq = MOD( kyr-1900. , 4. )
-      IF( zyq == 0.)   inc = 1.
+      !!                   ***  FUNCTION dayjul  ***
+      !!
+      !! Purpose :  compute the Julian day
+      !!----------------------------------------------------------------------
+      INTEGER,INTENT(in)    ::   kyr, kmonth, kday
+      !
+      INTEGER,DIMENSION(12) ::   idayt = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /)
+      INTEGER,DIMENSION(12) ::   idays
+      INTEGER               ::   inc, ji, zyq
+      REAL(wp)              ::   dayjul
+      !!----------------------------------------------------------------------
+      !
+      idays(1) = 0
+      idays(2) = 31
+      inc = 0.0_wp
+      zyq = MOD( kyr - 1900 , 4 )
+      IF( zyq == 0 ) inc = 1
       DO ji = 3, 12
-         idays(ji)=idayt(ji)+inc
+         idays(ji) = idayt(ji) + inc
       END DO
-      dayjul = idays(kmonth) + kday
+      dayjul = REAL( idays(kmonth) + kday, KIND=wp )
       !
    END FUNCTION dayjul
