New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3953 for branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90 – NEMO

Ignore:
Timestamp:
2013-07-03T13:41:32+02:00 (11 years ago)
Author:
gm
Message:

dev_r3858_NOC_ZTC, #863 : activate tide potential in filtered ssh case + style in tide modules

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90

    r3670 r3953  
    11MODULE tide_mod 
    2   !!================================================================================= 
    3   !!                       ***  MODULE  tide_mod  *** 
    4   !! Compute nodal modulations corrections and pulsations 
    5   !!================================================================================= 
    6   !!--------------------------------------------------------------------------------- 
    7   !!   OPA 9.0 , LODYC-IPSL  (2003) 
    8   !!--------------------------------------------------------------------------------- 
    9   USE dom_oce         ! ocean space and time domain 
    10   USE phycst 
    11   USE daymod 
    12  
    13   IMPLICIT NONE 
    14   PRIVATE 
    15  
    16   REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1, & 
    17        sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R, & 
    18        sh_I, sh_x1ra, sh_N 
    19  
    20   INTEGER,PUBLIC, PARAMETER ::   & 
    21        jpmax_harmo = 19             ! maximum number of harmonic 
    22  
    23   TYPE, PUBLIC ::    tide 
    24      CHARACTER(LEN=4)  :: cname_tide 
    25      REAL(wp) :: equitide 
    26      INTEGER  :: nutide 
    27      INTEGER  ::  nt,ns,nh,np,np1,shift 
    28      INTEGER  ::  nksi,nnu0,nnu1,nnu2,R 
    29      INTEGER  :: nformula 
    30   END TYPE tide 
    31  
    32   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave 
    33  
    34   !! * Accessibility 
    35   PUBLIC tide_harmo 
    36   PUBLIC nodal_factort 
    37   PUBLIC tide_init_Wave 
    38  
     2   !!====================================================================== 
     3   !!                       ***  MODULE  tide_mod  *** 
     4   !! Compute nodal modulations corrections and pulsations 
     5   !!====================================================================== 
     6   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code 
     7   !!---------------------------------------------------------------------- 
     8   USE dom_oce        ! ocean space and time domain 
     9   USE phycst         ! physical constant 
     10   USE daymod         ! calendar 
     11 
     12   IMPLICIT NONE 
     13   PRIVATE 
     14 
     15   PUBLIC   tide_harmo       ! called by tideini and diaharm modules 
     16   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules 
     17 
     18   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic 
     19 
     20   TYPE, PUBLIC ::    tide 
     21      CHARACTER(LEN=4) ::   cname_tide 
     22      REAL(wp)         ::   equitide 
     23      INTEGER          ::   nutide 
     24      INTEGER          ::   nt, ns, nh, np, np1, shift 
     25      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R 
     26      INTEGER          ::   nformula 
     27   END TYPE tide 
     28 
     29   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) ::   Wave   !: 
     30 
     31   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles 
     32   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   ! 
     33   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       ! 
     34 
     35   !!---------------------------------------------------------------------- 
     36   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)  
     37   !! $Id:$  
     38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     39   !!---------------------------------------------------------------------- 
    3940CONTAINS 
    4041 
    41   SUBROUTINE tide_init_Wave 
    42  
    43 #  include "tide.h90" 
    44  
    45   END SUBROUTINE tide_init_Wave 
    46  
    47   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc) 
    48  
    49     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    50          ktide      ! Indice of tidal constituents 
    51  
    52     INTEGER, INTENT( in ) :: & 
    53          kc         ! Total number of tidal constituents 
    54  
    55     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    56          pomega      ! pulsation in radians/s 
    57  
    58     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    59          pvt, &      ! 
    60          put, &      ! 
    61          pcor         ! 
    62  
    63     CALL astronomic_angle 
    64     CALL tide_pulse(pomega, ktide ,kc) 
    65     CALL tide_vuf( pvt, put, pcor, ktide ,kc) 
    66  
    67   END SUBROUTINE tide_harmo 
    68  
    69   SUBROUTINE astronomic_angle 
    70  
    71     !!---------------------------------------------------------------------- 
    72     !! 
    73     !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian 
    74     !!  century (e.g. time in days divide by 36525) 
    75     !!---------------------------------------------------------------------- 
    76  
    77     REAL(wp) ::  cosI,p,q,t2,t4,sin2I,s2,tgI2,P1,sh_tgn2,at1,at2 
    78     REAL(wp) :: zqy,zsy,zday,zdj,zhfrac 
    79  
    80     zqy=AINT((nyear-1901.)/4.) 
    81     zsy=nyear-1900. 
    82  
    83     zdj=dayjul(nyear,nmonth,nday) 
    84     zday=zdj+zqy-1. 
    85  
    86     zhfrac=nsec_day/3600. 
    87  
    88     !---------------------------------------------------------------------- 
    89     !  Sh_n Longitude of ascending lunar node 
    90     !---------------------------------------------------------------------- 
    91  
    92     sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad 
    93     !---------------------------------------------------------------------- 
    94     ! T mean solar angle (Greenwhich time) 
    95     !---------------------------------------------------------------------- 
    96     sh_T=(180.+zhfrac*(360./24.))*rad 
    97     !---------------------------------------------------------------------- 
    98     ! h mean solar Longitude 
    99     !---------------------------------------------------------------------- 
    100  
    101     sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad 
    102     !---------------------------------------------------------------------- 
    103     ! s mean lunar Longitude 
    104     !---------------------------------------------------------------------- 
    105  
    106     sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad 
    107     !---------------------------------------------------------------------- 
    108     ! p1 Longitude of solar perigee 
    109     !---------------------------------------------------------------------- 
    110  
    111     sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad 
    112     !---------------------------------------------------------------------- 
    113     ! p Longitude of lunar perigee 
    114     !---------------------------------------------------------------------- 
    115  
    116     sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 
    117  
    118     sh_N =mod(sh_N ,2*rpi) 
    119     sh_s =mod(sh_s ,2*rpi) 
    120     sh_h =mod(sh_h, 2*rpi) 
    121     sh_p =mod(sh_p, 2*rpi) 
    122     sh_p1=mod(sh_p1,2*rpi) 
    123  
    124     cosI=0.913694997 -0.035692561 *cos(sh_N) 
    125  
    126     sh_I=acos(cosI) 
    127  
    128     sin2I=sin(sh_I) 
    129     sh_tgn2=tan(sh_N/2.0) 
    130  
    131     at1=atan(1.01883*sh_tgn2) 
    132     at2=atan(0.64412*sh_tgn2) 
    133  
    134     sh_xi=-at1-at2+sh_N 
    135  
    136     if (sh_N > rpi) sh_xi=sh_xi-2.0*rpi 
    137  
    138     sh_nu=at1-at2 
    139  
    140     !---------------------------------------------------------------------- 
    141     ! For constituents l2 k1 k2 
    142     !---------------------------------------------------------------------- 
    143  
    144     tgI2=tan(sh_I/2.0) 
    145     P1=sh_p-sh_xi 
    146  
    147     t2=tgI2*tgI2 
    148     t4=t2*t2 
    149     sh_x1ra=sqrt(1.0-12.0*t2*cos(2.0*P1)+36.0*t4) 
    150  
    151     p=sin(2.0*P1) 
    152     q=1.0/(6.0*t2)-cos(2.0*P1) 
    153     sh_R=atan(p/q) 
    154  
    155     p=sin(2.0*sh_I)*sin(sh_nu) 
    156     q=sin(2.0*sh_I)*cos(sh_nu)+0.3347 
    157     sh_nuprim=atan(p/q) 
    158  
    159     s2=sin(sh_I)*sin(sh_I) 
    160     p=s2*sin(2.0*sh_nu) 
    161     q=s2*cos(2.0*sh_nu)+0.0727 
    162     sh_nusec=0.5*atan(p/q) 
    163  
    164   END SUBROUTINE astronomic_angle 
    165  
    166   SUBROUTINE tide_pulse( pomega, ktide ,kc) 
    167     !!---------------------------------------------------------------------- 
    168     !!                     ***  ROUTINE tide_pulse  *** 
    169     !!                       
    170     !! ** Purpose : Compute tidal frequencies 
    171     !! 
    172     !!---------------------------------------------------------------------- 
    173     !! * Arguments 
    174     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    175          ktide      ! Indice of tidal constituents 
    176  
    177     INTEGER, INTENT( in ) :: & 
    178          kc         ! Total number of tidal constituents 
    179  
    180     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    181          pomega      ! pulsation in radians/s 
    182  
    183     !! * Local declarations 
    184     INTEGER :: jh 
    185     REAL(wp) :: zscale  =  36525*24.0 
    186     REAL(wp) :: zomega_T=  13149000.0 
    187     REAL(wp) :: zomega_s=    481267.892 
    188     REAL(wp) :: zomega_h=     36000.76892 
    189     REAL(wp) :: zomega_p=      4069.0322056 
    190     REAL(wp) :: zomega_n=      1934.1423972 
    191     REAL(wp) :: zomega_p1=        1.719175 
    192     !!---------------------------------------------------------------------- 
    193  
    194     DO jh=1,kc 
    195        pomega(jh) = zomega_T * Wave(ktide(jh))%nT & 
    196             + zomega_s * Wave(ktide(jh))%ns & 
    197             + zomega_h * Wave(ktide(jh))%nh & 
    198             + zomega_p * Wave(ktide(jh))%np & 
    199             + zomega_p1* Wave(ktide(jh))%np1 
    200        pomega(jh) = (pomega(jh)/zscale)*rad/3600. 
    201     END DO 
    202  
    203   END SUBROUTINE tide_pulse 
    204  
    205   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc) 
    206     !!---------------------------------------------------------------------- 
    207     !!                     ***  ROUTINE tide_vuf  *** 
    208     !!                       
    209     !! ** Purpose : Compute nodal modulation corrections 
    210     !! 
    211     !! ** Outputs : 
    212     !!          vt: Pase of tidal potential relative to Greenwich (radians) 
    213     !!          ut: Phase correction u due to nodal motion (radians) 
    214     !!          ft: Nodal correction factor 
    215     !! 
    216     !! ** Inputs : 
    217     !!          tname: array of constituents names (dimension<=nc)  
    218     !!             nc: number of constituents 
    219     !!    
    220     !!---------------------------------------------------------------------- 
    221     !! * Arguments 
    222     INTEGER, DIMENSION(kc), INTENT( in ) ::   & 
    223          ktide      ! Indice of tidal constituents 
    224     INTEGER, INTENT( in ) :: & 
    225          kc         ! Total number of tidal constituents 
    226     REAL (wp), DIMENSION(kc), INTENT( out ) ::   & 
    227          pvt, &      ! 
    228          put, &      ! 
    229          pcor         ! 
    230     !! * Local declarations 
    231     INTEGER :: jh 
    232     !!---------------------------------------------------------------------- 
    233  
    234     DO jh =1,kc 
    235        !  Phase of the tidal potential relative to the Greenwhich  
    236        !  meridian (e.g. the position of the fictuous celestial body). Units are 
    237        !  radian: 
    238        pvt(jh) = sh_T *Wave(ktide(jh))%nT        & 
    239             +sh_s *Wave(ktide(jh))%ns        & 
    240             +sh_h *Wave(ktide(jh))%nh        & 
    241             +sh_p *Wave(ktide(jh))%np        & 
    242             +sh_p1*Wave(ktide(jh))%np1       & 
    243             +Wave(ktide(jh))%shift*rad 
     42   SUBROUTINE tide_init_Wave 
     43#     include "tide.h90" 
     44   END SUBROUTINE tide_init_Wave 
     45 
     46 
     47   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc) 
     48      !!---------------------------------------------------------------------- 
     49      !!---------------------------------------------------------------------- 
     50      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents 
     51      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents 
     52      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega           ! pulsation in radians/s 
     53      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   ! 
     54      !!---------------------------------------------------------------------- 
     55      ! 
     56      CALL astronomic_angle 
     57      CALL tide_pulse( pomega, ktide ,kc ) 
     58      CALL tide_vuf  ( pvt, put, pcor, ktide ,kc ) 
     59      ! 
     60   END SUBROUTINE tide_harmo 
     61 
     62 
     63   SUBROUTINE astronomic_angle 
     64      !!---------------------------------------------------------------------- 
     65      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian 
     66      !!  century (e.g. time in days divide by 36525) 
     67      !!---------------------------------------------------------------------- 
     68      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2 
     69      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac 
     70      !!---------------------------------------------------------------------- 
     71      ! 
     72      zqy = AINT( (nyear-1901.)/4. ) 
     73      zsy = nyear - 1900. 
     74      ! 
     75      zdj  = dayjul( nyear, nmonth, nday ) 
     76      zday = zdj + zqy - 1. 
     77      ! 
     78      zhfrac = nsec_day / 3600. 
     79      ! 
     80      !---------------------------------------------------------------------- 
     81      !  Sh_n Longitude of ascending lunar node 
     82      !---------------------------------------------------------------------- 
     83      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad 
     84      !---------------------------------------------------------------------- 
     85      ! T mean solar angle (Greenwhich time) 
     86      !---------------------------------------------------------------------- 
     87      sh_T=(180.+zhfrac*(360./24.))*rad 
     88      !---------------------------------------------------------------------- 
     89      ! h mean solar Longitude 
     90      !---------------------------------------------------------------------- 
     91      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad 
     92      !---------------------------------------------------------------------- 
     93      ! s mean lunar Longitude 
     94      !---------------------------------------------------------------------- 
     95      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad 
     96      !---------------------------------------------------------------------- 
     97      ! p1 Longitude of solar perigee 
     98      !---------------------------------------------------------------------- 
     99      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad 
     100      !---------------------------------------------------------------------- 
     101      ! p Longitude of lunar perigee 
     102      !---------------------------------------------------------------------- 
     103      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 
     104 
     105      sh_N = MOD( sh_N ,2*rpi ) 
     106      sh_s = MOD( sh_s ,2*rpi ) 
     107      sh_h = MOD( sh_h, 2*rpi ) 
     108      sh_p = MOD( sh_p, 2*rpi ) 
     109      sh_p1= MOD( sh_p1,2*rpi ) 
     110 
     111      cosI = 0.913694997 -0.035692561 *cos(sh_N) 
     112 
     113      sh_I = ACOS( cosI ) 
     114 
     115      sin2I   = sin(sh_I) 
     116      sh_tgn2 = tan(sh_N/2.0) 
     117 
     118      at1=atan(1.01883*sh_tgn2) 
     119      at2=atan(0.64412*sh_tgn2) 
     120 
     121      sh_xi=-at1-at2+sh_N 
     122 
     123      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi 
     124 
     125      sh_nu = at1 - at2 
     126 
     127      !---------------------------------------------------------------------- 
     128      ! For constituents l2 k1 k2 
     129      !---------------------------------------------------------------------- 
     130 
     131      tgI2 = tan(sh_I/2.0) 
     132      P1   = sh_p-sh_xi 
     133 
     134      t2 = tgI2*tgI2 
     135      t4 = t2*t2 
     136      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 ) 
     137 
     138      p = sin(2.0*P1) 
     139      q = 1.0/(6.0*t2)-cos(2.0*P1) 
     140      sh_R = atan(p/q) 
     141 
     142      p = sin(2.0*sh_I)*sin(sh_nu) 
     143      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347 
     144      sh_nuprim = atan(p/q) 
     145 
     146      s2 = sin(sh_I)*sin(sh_I) 
     147      p  = s2*sin(2.0*sh_nu) 
     148      q  = s2*cos(2.0*sh_nu)+0.0727 
     149      sh_nusec = 0.5*atan(p/q) 
     150      ! 
     151   END SUBROUTINE astronomic_angle 
     152 
     153 
     154   SUBROUTINE tide_pulse( pomega, ktide ,kc ) 
     155      !!---------------------------------------------------------------------- 
     156      !!                     ***  ROUTINE tide_pulse  *** 
     157      !!                       
     158      !! ** Purpose : Compute tidal frequencies 
     159      !!---------------------------------------------------------------------- 
     160      INTEGER                , INTENT(in ) ::   kc       ! Total number of tidal constituents 
     161      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide    ! Indice of tidal constituents 
     162      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega   ! pulsation in radians/s 
     163      ! 
     164      INTEGER  ::   jh 
     165      REAL(wp) ::   zscale 
     166      REAL(wp) ::   zomega_T =  13149000.0_wp 
     167      REAL(wp) ::   zomega_s =    481267.892_wp 
     168      REAL(wp) ::   zomega_h =     36000.76892_wp 
     169      REAL(wp) ::   zomega_p =      4069.0322056_wp 
     170      REAL(wp) ::   zomega_n =      1934.1423972_wp 
     171      REAL(wp) ::   zomega_p1=         1.719175_wp 
     172      !!---------------------------------------------------------------------- 
     173      ! 
     174      zscale =  rad / ( 36525._wp * 86400._wp )  
     175      ! 
     176      DO jh = 1, kc 
     177         pomega(jh) = (  zomega_T * Wave( ktide(jh) )%nT   & 
     178            &          + zomega_s * Wave( ktide(jh) )%ns   & 
     179            &          + zomega_h * Wave( ktide(jh) )%nh   & 
     180            &          + zomega_p * Wave( ktide(jh) )%np   & 
     181            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale 
     182      END DO 
     183      ! 
     184   END SUBROUTINE tide_pulse 
     185 
     186 
     187   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc ) 
     188      !!---------------------------------------------------------------------- 
     189      !!                     ***  ROUTINE tide_vuf  *** 
     190      !!                       
     191      !! ** Purpose : Compute nodal modulation corrections 
     192      !! 
     193      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians) 
     194      !!              ut: Phase correction u due to nodal motion (radians) 
     195      !!              ft: Nodal correction factor 
     196      !!---------------------------------------------------------------------- 
     197      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents 
     198      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents 
     199      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   ! 
     200      ! 
     201      INTEGER ::   jh   ! dummy loop index 
     202      !!---------------------------------------------------------------------- 
     203      ! 
     204      DO jh = 1, kc 
     205         !  Phase of the tidal potential relative to the Greenwhich  
     206         !  meridian (e.g. the position of the fictuous celestial body). Units are radian: 
     207         pvt(jh) = sh_T * Wave( ktide(jh) )%nT    & 
     208            &    + sh_s * Wave( ktide(jh) )%ns    & 
     209            &    + sh_h * Wave( ktide(jh) )%nh    & 
     210            &    + sh_p * Wave( ktide(jh) )%np    & 
     211            &    + sh_p1* Wave( ktide(jh) )%np1   & 
     212            &    +        Wave( ktide(jh) )%shift * rad 
     213         ! 
     214         !  Phase correction u due to nodal motion. Units are radian: 
     215         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   & 
     216            &    + sh_nu     * Wave( ktide(jh) )%nnu0   & 
     217            &    + sh_nuprim * Wave( ktide(jh) )%nnu1   & 
     218            &    + sh_nusec  * Wave( ktide(jh) )%nnu2   & 
     219            &    + sh_R      * Wave( ktide(jh) )%R 
     220 
     221         !  Nodal correction factor: 
     222         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula ) 
     223      END DO 
     224      ! 
     225   END SUBROUTINE tide_vuf 
     226 
     227 
     228   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf ) 
     229      !!---------------------------------------------------------------------- 
     230      !!---------------------------------------------------------------------- 
     231      INTEGER, INTENT(in) :: kformula 
     232      ! 
     233      REAL(wp) :: zf 
     234      REAL(wp) :: zs, zf1, zf2 
     235      !!---------------------------------------------------------------------- 
     236      ! 
     237      SELECT CASE( kformula ) 
     238      ! 
     239      CASE( 0 )                  !==  formule 0, solar waves 
     240         zf = 1.0 
     241         ! 
     242      CASE( 1 )                  !==  formule 1, compound waves (78 x 78) 
     243         zf=nodal_factort(78) 
     244         zf = zf * zf 
     245         ! 
     246      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78)  
     247       zf1= nodal_factort(78) 
     248       zf = nodal_factort( 0) 
     249       zf = zf1 * zf 
    244250       ! 
    245        !  Phase correction u due to nodal motion. Units are radian: 
    246        put(jh) = sh_xi    *Wave(ktide(jh))%nksi  & 
    247             +sh_nu    *Wave(ktide(jh))%nnu0  & 
    248             +sh_nuprim*Wave(ktide(jh))%nnu1  & 
    249             +sh_nusec *Wave(ktide(jh))%nnu2  & 
    250             +sh_R     *Wave(ktide(jh))%R 
    251  
    252        !  Nodal correction factor: 
    253        pcor(jh) = nodal_factort(Wave(ktide(jh))%nformula) 
    254     END DO 
    255  
    256   END SUBROUTINE tide_vuf 
    257  
    258   recursive function nodal_factort(kformula) result (zf) 
    259     !!---------------------------------------------------------------------- 
    260     INTEGER, INTENT(IN) :: kformula 
    261     REAL(wp) :: zf 
    262     REAL(wp) :: zs,zf1,zf2 
    263  
    264     SELECT CASE (kformula) 
    265  
    266        !!  formule 0, solar waves 
    267  
    268     case ( 0 ) 
    269        zf=1.0 
    270  
    271        !! formule 1, compound waves (78 x 78) 
    272  
    273     case ( 1 ) 
    274        zf=nodal_factort(78) 
    275        zf=zf*zf 
    276  
    277        !! formule 2, compound waves (78 x 0)  ===  (78)  
    278  
    279     case ( 2 ) 
    280        zf1=nodal_factort(78) 
    281        zf=nodal_factort(0) 
    282        zf=zf1*zf 
    283  
    284        !! formule 4,  compound waves (78 x 235)  
    285  
    286     case ( 4 ) 
    287        zf1=nodal_factort(78) 
    288        zf=nodal_factort(235) 
    289        zf=zf1*zf 
    290  
    291        !! formule 5,  compound waves (78 *78 x 235) 
    292  
    293     case ( 5 ) 
    294        zf1=nodal_factort(78) 
    295        zf=nodal_factort(235) 
    296        zf=zf*zf1*zf1 
    297  
    298        !! formule 6,  compound waves (78 *78 x 0) 
    299  
    300     case ( 6 ) 
    301        zf1=nodal_factort(78) 
    302        zf=nodal_factort(0) 
    303        zf=zf*zf1*zf1  
    304  
    305        !! formule 7, compound waves (75 x 75) 
    306  
    307     case ( 7 ) 
    308        zf=nodal_factort(75) 
    309        zf=zf*zf 
    310  
    311        !! formule 8,  compound waves (78 x 0 x 235) 
    312  
    313     case ( 8 ) 
    314        zf=nodal_factort(78) 
    315        zf1=nodal_factort(0) 
    316        zf2=nodal_factort(235) 
    317        zf=zf*zf1*zf2 
    318  
    319        !! formule 9,  compound waves (78 x 0 x 227) 
    320  
    321     case ( 9 ) 
    322        zf=nodal_factort(78) 
    323        zf1=nodal_factort(0) 
    324        zf2=nodal_factort(227) 
    325        zf=zf*zf1*zf2 
    326  
    327        !! formule 10,  compound waves (78 x 227) 
    328  
    329     case ( 10 ) 
    330        zf=nodal_factort(78) 
    331        zf1=nodal_factort(227) 
    332        zf=zf*zf1 
    333  
    334        !! formule 11,  compound waves (75 x 0) 
    335  
    336     case ( 11 ) 
    337        zf=nodal_factort(75) 
    338        zf=nodal_factort(0) 
    339        zf=zf*zf1 
    340  
    341        !! formule 12,  compound waves (78 x 78 x 78 x 0)  
    342  
    343     case ( 12 ) 
    344        zf1=nodal_factort(78) 
    345        zf=nodal_factort(0) 
    346        zf=zf*zf1*zf1*zf1 
    347  
    348        !! formule 13, compound waves (78 x 75) 
    349  
    350     case ( 13 ) 
    351        zf1=nodal_factort(78) 
    352        zf=nodal_factort(75) 
    353        zf=zf*zf1 
    354  
    355        !! formule 14, compound waves (235 x 0)  ===  (235) 
    356  
    357     case ( 14 ) 
    358        zf=nodal_factort(235) 
    359        zf1=nodal_factort(0) 
    360        zf=zf*zf1 
    361  
    362        !! formule 15, compound waves (235 x 75)  
    363  
    364     case ( 15 ) 
    365        zf=nodal_factort(235) 
    366        zf1=nodal_factort(75) 
    367        zf=zf*zf1 
    368  
    369        !! formule 16, compound waves (78 x 0 x 0)  ===  (78) 
    370  
    371     case ( 16 ) 
    372        zf=nodal_factort(78) 
    373        zf1=nodal_factort(0) 
    374        zf=zf*zf1*zf1 
    375  
    376        !! formule 17,  compound waves (227 x 0)  
    377  
    378     case ( 17 ) 
    379        zf1=nodal_factort(227) 
    380        zf=nodal_factort(0) 
    381        zf=zf*zf1 
    382  
    383        !! formule 18,  compound waves (78 x 78 x 78 ) 
    384  
    385     case ( 18 )  
    386        zf1=nodal_factort(78) 
    387        zf=zf1*zf1*zf1 
    388  
    389        !! formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78) 
    390  
    391     case ( 19 ) 
    392        zf=nodal_factort(78) 
    393        zf1=nodal_factort(0) 
    394        zf=zf*zf1*zf1 
    395  
    396        !! formule 73 
    397  
    398     case ( 73 ) 
    399        zs=sin(sh_I) 
    400        zf=(2./3.-zs*zs)/0.5021 
    401  
    402        !! formule 74 
    403  
    404     case ( 74 ) 
    405        zs=sin(sh_I) 
    406        zf=zs*zs/0.1578 
    407  
    408        !! formule 75 
    409  
    410     case ( 75 ) 
    411        zs=cos (sh_I/2) 
    412        zf=sin (sh_I)*zs*zs/0.3800 
    413  
    414        !! formule 76 
    415  
    416     case ( 76 ) 
    417        zf=sin (2*sh_I)/0.7214 
    418  
    419        !! formule 77 
    420  
    421     case ( 77 ) 
    422        zs=sin (sh_I/2) 
    423        zf=sin (sh_I)*zs*zs/0.0164 
    424  
    425        !! formule 78 
    426  
    427     case ( 78 ) 
    428        zs=cos (sh_I/2) 
    429        zf=zs*zs*zs*zs/0.9154 
    430  
    431        !! formule 79 
    432  
    433     case ( 79 ) 
    434        zs=sin(sh_I) 
    435        zf=zs*zs/0.1565 
    436  
    437        !! formule 144 
    438  
    439     case ( 144 ) 
    440        zs=sin (sh_I/2) 
    441        zf=(1-10*zs*zs+15*zs*zs*zs*zs)*cos(sh_I/2)/0.5873 
    442  
    443        !! formule 149 
    444  
    445     case ( 149 ) 
    446        zs=cos (sh_I/2) 
    447        zf=zs*zs*zs*zs*zs*zs/0.8758 
    448  
    449        !! formule 215 
    450  
    451     case ( 215 ) 
    452        zs=cos (sh_I/2) 
    453        zf=zs*zs*zs*zs/0.9154*sh_x1ra 
    454  
    455        !! formule 227  
    456  
    457     case ( 227 ) 
    458        zs=sin (2*sh_I) 
    459        zf=sqrt (0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006) 
    460  
    461        !! formule 235  
    462  
    463     case ( 235 ) 
    464        zs=sin (sh_I) 
    465        zf=sqrt (19.0444*zs*zs*zs*zs+2.7702*zs*zs*cos (2*sh_nu)+.0981) 
    466  
    467     END SELECT 
    468  
    469   end function nodal_factort 
    470  
    471   function dayjul(kyr,kmonth,kday) 
    472     ! 
    473     !*** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE) 
    474     ! 
    475     INTEGER,INTENT(IN) :: kyr,kmonth,kday 
    476     INTEGER,DIMENSION(12) ::  idayt,idays 
    477     INTEGER :: inc,ji 
    478     REAL(wp) :: dayjul,zyq 
    479  
    480     DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 
    481     idays(1)=0. 
    482     idays(2)=31. 
    483     inc=0. 
    484     zyq=MOD((kyr-1900.),4.) 
    485     IF(zyq .eq. 0.) inc=1. 
    486     DO ji=3,12 
    487        idays(ji)=idayt(ji)+inc 
    488     END DO 
    489     dayjul=idays(kmonth)+kday 
    490  
    491   END FUNCTION dayjul 
    492  
     251      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235)  
     252         zf1 = nodal_factort( 78) 
     253         zf  = nodal_factort(235) 
     254         zf  = zf1 * zf 
     255         ! 
     256      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235) 
     257         zf1 = nodal_factort( 78) 
     258         zf  = nodal_factort(235) 
     259         zf  = zf * zf1 * zf1 
     260         ! 
     261      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0) 
     262         zf1 = nodal_factort(78) 
     263         zf  = nodal_factort( 0) 
     264         zf  = zf * zf1 * zf1  
     265         ! 
     266      CASE( 7 )                  !==  formule 7, compound waves (75 x 75) 
     267         zf = nodal_factort(75) 
     268         zf = zf * zf 
     269         ! 
     270      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235) 
     271         zf  = nodal_factort( 78) 
     272         zf1 = nodal_factort(  0) 
     273         zf2 = nodal_factort(235) 
     274         zf  = zf * zf1 * zf2 
     275         ! 
     276      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227) 
     277         zf  = nodal_factort( 78) 
     278         zf1 = nodal_factort(  0) 
     279         zf2 = nodal_factort(227) 
     280         zf  = zf * zf1 * zf2 
     281         ! 
     282      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227) 
     283         zf  = nodal_factort( 78) 
     284         zf1 = nodal_factort(227) 
     285         zf  = zf * zf1 
     286         ! 
     287      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0) 
     288!!gm bug???? zf 2 fois ! 
     289         zf = nodal_factort(75) 
     290         zf = nodal_factort( 0) 
     291         zf = zf * zf1 
     292         ! 
     293      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0)  
     294         zf1 = nodal_factort(78) 
     295         zf  = nodal_factort( 0) 
     296         zf  = zf * zf1 * zf1 * zf1 
     297         ! 
     298      CASE( 13 )                 !==  formule 13, compound waves (78 x 75) 
     299         zf1 = nodal_factort(78) 
     300         zf  = nodal_factort(75) 
     301         zf  = zf * zf1 
     302         ! 
     303      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235) 
     304         zf  = nodal_factort(235) 
     305         zf1 = nodal_factort(  0) 
     306         zf  = zf * zf1 
     307         ! 
     308      CASE( 15 )                 !==  formule 15, compound waves (235 x 75)  
     309         zf  = nodal_factort(235) 
     310         zf1 = nodal_factort( 75) 
     311         zf  = zf * zf1 
     312         ! 
     313      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78) 
     314         zf  = nodal_factort(78) 
     315         zf1 = nodal_factort( 0) 
     316         zf  = zf * zf1 * zf1 
     317         ! 
     318      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0)  
     319         zf1 = nodal_factort(227) 
     320         zf  = nodal_factort(  0) 
     321         zf  = zf * zf1 
     322         ! 
     323      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 ) 
     324         zf1 = nodal_factort(78) 
     325         zf  = zf1 * zf1 * zf1 
     326         ! 
     327      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78) 
     328!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing 
     329         zf  = nodal_factort(78) 
     330         zf1 = nodal_factort( 0) 
     331         zf = zf * zf1 * zf1 
     332         ! 
     333      CASE( 73 )                 !==  formule 73 
     334         zs = sin(sh_I) 
     335         zf = (2./3.-zs*zs)/0.5021 
     336         ! 
     337      CASE( 74 )                 !==  formule 74 
     338         zs = sin(sh_I) 
     339         zf = zs * zs / 0.1578 
     340         ! 
     341      CASE( 75 )                 !==  formule 75 
     342         zs = cos(sh_I/2) 
     343         zf = sin(sh_I) * zs * zs / 0.3800 
     344         ! 
     345      CASE( 76 )                 !==  formule 76 
     346         zf = sin(2*sh_I) / 0.7214 
     347         ! 
     348      CASE( 77 )                 !==  formule 77 
     349         zs = sin(sh_I/2) 
     350         zf = sin(sh_I) * zs * zs / 0.0164 
     351         ! 
     352      CASE( 78 )                 !==  formule 78 
     353         zs = cos(sh_I/2) 
     354         zf = zs * zs * zs * zs / 0.9154 
     355         ! 
     356      CASE( 79 )                 !==  formule 79 
     357         zs = sin(sh_I) 
     358         zf = zs * zs / 0.1565 
     359         ! 
     360      CASE( 144 )                !==  formule 144 
     361         zs = sin(sh_I/2) 
     362         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873 
     363         ! 
     364      CASE( 149 )                !==  formule 149 
     365         zs = cos(sh_I/2) 
     366         zf = zs*zs*zs*zs*zs*zs / 0.8758 
     367         ! 
     368      CASE( 215 )                !==  formule 215 
     369         zs = cos(sh_I/2) 
     370         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra 
     371         ! 
     372      CASE( 227 )                !==  formule 227  
     373         zs = sin(2*sh_I) 
     374         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 ) 
     375         ! 
     376      CASE ( 235 )               !==  formule 235  
     377         zs = sin(sh_I) 
     378         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 ) 
     379         ! 
     380      END SELECT 
     381      ! 
     382   END FUNCTION nodal_factort 
     383 
     384 
     385   FUNCTION dayjul( kyr, kmonth, kday ) 
     386      !!---------------------------------------------------------------------- 
     387      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE) 
     388      !!---------------------------------------------------------------------- 
     389      INTEGER,INTENT(in) ::   kyr, kmonth, kday 
     390      ! 
     391      INTEGER,DIMENSION(12) ::  idayt, idays 
     392      INTEGER  ::   inc, ji 
     393      REAL(wp) ::   dayjul, zyq 
     394      ! 
     395      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 
     396      !!---------------------------------------------------------------------- 
     397      ! 
     398      idays(1) = 0. 
     399      idays(2) = 31. 
     400      inc = 0. 
     401      zyq = MOD( kyr-1900. , 4. ) 
     402      IF( zyq == 0.)   inc = 1. 
     403      DO ji = 3, 12 
     404         idays(ji)=idayt(ji)+inc 
     405      END DO 
     406      dayjul = idays(kmonth) + kday 
     407      ! 
     408   END FUNCTION dayjul 
     409 
     410   !!====================================================================== 
    493411END MODULE tide_mod 
Note: See TracChangeset for help on using the changeset viewer.