Changeset 12042 for NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90

Timestamp:

2019-12-03T16:41:39+01:00 (4 years ago)

Author:

smueller

Message:

Reduction of the number of real-valued literals used in subroutine astronomic_angle of module tide_mod (ticket #2194)

File:

• NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90 (modified) (2 diffs)

NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90

@@ -67 +67 @@
    REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
    REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !
+
+   ! Longitudes on 1 Jan 1900, 00h and angular velocities (units of deg and
+   ! deg/h, respectively. The values of these module variables have been copied
+   ! from subroutine astronomic_angle of the version of this module used in
+   ! version 4.0 of NEMO
+   REAL(wp) ::   rlon00_N  =  259.1560564_wp               ! Longitude of ascending lunar node
+   REAL(wp) ::   romega_N  = -.0022064139_wp
+   REAL(wp) ::   rlon00_T  =  180.0_wp                     ! Mean solar angle (GMT)
+   REAL(wp) ::   romega_T  =  15.0_wp
+   REAL(wp) ::   rlon00_h  =  280.1895014_wp               ! Mean solar Longitude
+   REAL(wp) ::   romega_h  =  .0410686387_wp
+   REAL(wp) ::   rlon00_s  =  277.0256206_wp               ! Mean lunar Longitude
+   REAL(wp) ::   romega_s  =  .549016532_wp
+   REAL(wp) ::   rlon00_p1 =  281.2208569_wp               ! Longitude of solar perigee
+   REAL(wp) ::   romega_p1 =  .000001961_wp
+   REAL(wp) ::   rlon00_p  =  334.3837214_wp               ! Longitude of lunar perigee
+   REAL(wp) ::   romega_p  =  .004641834_wp
 
    !!----------------------------------------------------------------------
@@ -423 +440 @@
    SUBROUTINE astronomic_angle(psec_day)
       !!----------------------------------------------------------------------
-      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
-      !!  century (e.g. time in days divide by 36525)
+      !!                 ***  ROUTINE astronomic_angle  ***
+      !!                      
+      !! ** Purpose : Compute astronomic angles
       !!----------------------------------------------------------------------
       INTEGER  ::   psec_day !   Number of seconds from midnight
       REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
-      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
-      !!----------------------------------------------------------------------
-      !
-      zqy = AINT( (nyear-1901.)/4. )
-      zsy = nyear - 1900.
+      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac, zt
+      !!----------------------------------------------------------------------
+      !
+      ! Computation of the time from 1 Jan 1900, 00h in years
+      zqy = AINT( (nyear - 1901.0_wp) / 4.0_wp )
+      zsy = nyear - 1900.0_wp
       !
       zdj  = dayjul( nyear, nmonth, nday )
-      zday = zdj + zqy - 1.
-      !
-      zhfrac = psec_day / 3600.
-      !
-      !----------------------------------------------------------------------
-      !  Sh_n Longitude of ascending lunar node
-      !----------------------------------------------------------------------
-      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
-      !----------------------------------------------------------------------
-      ! T mean solar angle (Greenwhich time)
-      !----------------------------------------------------------------------
-      sh_T=(180.+zhfrac*(360./24.))*rad
-      !----------------------------------------------------------------------
-      ! h mean solar Longitude
-      !----------------------------------------------------------------------
-      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
-      !----------------------------------------------------------------------
-      ! s mean lunar Longitude
-      !----------------------------------------------------------------------
-      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
-      !----------------------------------------------------------------------
-      ! p1 Longitude of solar perigee
-      !----------------------------------------------------------------------
-      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
-      !----------------------------------------------------------------------
-      ! p Longitude of lunar perigee
-      !----------------------------------------------------------------------
-      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
-
-      sh_N = MOD( sh_N ,2*rpi )
-      sh_s = MOD( sh_s ,2*rpi )
-      sh_h = MOD( sh_h, 2*rpi )
-      sh_p = MOD( sh_p, 2*rpi )
-      sh_p1= MOD( sh_p1,2*rpi )
+      zday = zdj + zqy - 1.0_wp
+      !
+      zhfrac = psec_day / 3600.0_wp
+      !
+      zt = zsy * 365.0_wp * 24.0_wp + zday * 24.0_wp + zhfrac
+      !
+      ! Longitude of ascending lunar node
+      sh_N  = ( rlon00_N  + romega_N * zt     ) * rad
+      sh_N = MOD( sh_N,  2*rpi )
+      ! Mean solar angle (Greenwhich time)
+      sh_T  = ( rlon00_T  + romega_T * zhfrac ) * rad
+      ! Mean solar Longitude
+      sh_h  = ( rlon00_h  + romega_h * zt     ) * rad
+      sh_h = MOD( sh_h,  2*rpi )
+      ! Mean lunar Longitude
+      sh_s  = ( rlon00_s  + romega_s * zt     ) * rad
+      sh_s = MOD( sh_s,  2*rpi )
+      ! Longitude of solar perigee
+      sh_p1 = ( rlon00_p1 + romega_p1 * zt    ) * rad
+      sh_p1= MOD( sh_p1, 2*rpi )
+      ! Longitude of lunar perigee
+      sh_p  = ( rlon00_p  + romega_p  * zt    ) * rad
+      sh_p = MOD( sh_p,  2*rpi )
 
       cosI = 0.913694997 -0.035692561 *cos(sh_N)