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 15490 – NEMO

Changeset 15490


Ignore:
Timestamp:
2021-11-10T12:33:18+01:00 (6 months ago)
Author:
hadjt
Message:

SBC/tide_mod.F90

Adding updated tide module to handle 360 day tides. This is working, but needs some more work.

Care need for compress option, as unpublished, work in progress. Worked needed to set phase offsets (in v0tide).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC/tide_mod.F90

    r15456 r15490  
    66   !! History :  1.0  !  2007  (O. Le Galloudec)  Original code 
    77   !!---------------------------------------------------------------------- 
     8   !JT USE oce             ! ocean dynamics and tracers 
     9   USE lib_mpp         ! distributed memory computing library 
     10 
    811   USE dom_oce        ! ocean space and time domain 
    912   USE phycst         ! physical constant 
    1013   USE daymod         ! calendar 
     14   USE in_out_manager ! I/O units 
     15   USE ioipsl  , ONLY :   ymds2ju      ! for calendar 
     16 
    1117 
    1218   IMPLICIT NONE 
     
    3440   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       ! 
    3541 
     42   !!JT 
     43   INTEGER(KIND=8)  ::  days_since_origin 
     44   LOGICAL  ::   ln_astro_verbose  
     45   !LOGICAL  ::   ln_tide_360_cal  
     46   !LOGICAL  ::   ln_tide_drift_time_cont_manual 
     47   LOGICAL  ::   ln_tide_drift                  ! Do we want to run with "drifting" tides? (Namelist) 
     48   LOGICAL  ::   ln_tide_compress               ! Do we want to run with "compressed" tides? (Namelist) 
     49   INTEGER  ::   nn_tide_orig_yr,nn_tide_orig_mn,nn_tide_orig_dy            !JT 
     50 
     51   !!JT 
     52 
    3653   !!---------------------------------------------------------------------- 
    3754   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4764 
    4865   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc) 
     66 
     67      !! Externally called by sbctide.F90/sbc_tide 
     68      !! Externally named: omega_tide, v0tide, utide, ftide, ntide, nb_harmo 
    4969      !!---------------------------------------------------------------------- 
    5070      !!---------------------------------------------------------------------- 
     
    5575      !!---------------------------------------------------------------------- 
    5676      ! 
     77 
     78      INTEGER                              ::   ios 
     79 
     80 
     81      ln_tide_drift = .FALSE. 
     82      ln_tide_compress = .FALSE. 
     83 
     84      NAMELIST/nam_tides360/ ln_tide_drift,ln_tide_compress,ln_astro_verbose,& 
     85        & nn_tide_orig_yr,nn_tide_orig_mn,nn_tide_orig_dy 
     86 
     87      ! read in Namelist.  
     88      !!---------------------------------------------------------------------- 
     89      ! 
     90      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics 
     91      READ   ( numnam_ref, nam_tides360, IOSTAT=ios, ERR= 901 ) 
     92901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tides360 in reference namelist' ) 
     93 
     94      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics 
     95      READ  ( numnam_cfg, nam_tides360, IOSTAT = ios, ERR = 902 ) 
     96902   IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tides360 in configuration namelist' ) 
     97      IF(lwm) WRITE ( numond, nam_tides360 ) 
     98 
     99 
     100      IF( lwp ) THEN 
     101        WRITE(numout,*) " " 
     102        WRITE(numout,*) "tide_harmo: nam_tides360 - 360 day tides " 
     103        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
     104        WRITE(numout,*) "       tides360: allow tides to drift through year: ln_tide_drift = ",ln_tide_drift 
     105        WRITE(numout,*) "       tides360: Compress tides, so around a 360 day year: ln_tide_compress = ",ln_tide_compress 
     106        WRITE(numout,*) "       tides360:           USE ln_tide_compress  WITH CARE. INCOMPLETE." 
     107        WRITE(numout,*) "       tides360: Increase output verbosity: ln_astro_verbose = ",ln_astro_verbose 
     108        !WRITE(numout,*) "       tides360: Calculate time between origin and gregorian and 360 manually: ln_tide_drift_time_cont_manual = ",ln_tide_drift_time_cont_manual 
     109        WRITE(numout,*) "       tides360: 360 day origin date year: nn_tide_orig_yr = ",nn_tide_orig_yr 
     110        WRITE(numout,*) "       tides360: 360 day origin date month: nn_tide_orig_mn = ",nn_tide_orig_mn 
     111        WRITE(numout,*) "       tides360: 360 day origin date day: nn_tide_orig_dy = ",nn_tide_orig_dy 
     112        WRITE(numout,*) " " 
     113      ENDIF 
     114 
     115  
     116      IF( nleapy == 30 ) THEN 
     117          IF ( ln_tide_drift .AND. ln_tide_compress ) THEN 
     118              CALL ctl_stop( 'tide_harmo: nam_tides360: if 360 day calendar ln_tide_drift and ln_tide_compress cannot be true' ) 
     119          ENDIF 
     120           
     121 
     122          IF ( ln_tide_drift   ) THEN 
     123              WRITE(numout,*) "       tides360: Tides continuous so equinoctal tides drift through the year," 
     124              WRITE(numout,*) "                 as the S2-K2 beating occurs 5 days later every year." 
     125          ENDIF 
     126 
     127          IF ( ln_tide_compress   ) THEN 
     128              WRITE(numout,*) "       tides360: The Tropical Year (and so some tidal periods) are compressed," 
     129              WRITE(numout,*) "                 so the tides repeat with an annual cycle, so the " 
     130              WRITE(numout,*) "                 the S2-K2 beating is fixed relative to the calendar, but the " 
     131              WRITE(numout,*) "                 M2 period varies slightly." 
     132              WRITE(numout,*) "                 Use with care, as this requires more work." 
     133          ENDIF 
     134 
     135          IF ( ( .NOT. ln_tide_drift  ) .AND. ( .NOT. ln_tide_compress ) ) THEN 
     136              WRITE(numout,*) "       tides360: Use the default NEMO tide code, where the tides are reset " 
     137              WRITE(numout,*) "                 at the beginning of each month, leading to a slight discontinuity" 
     138              WRITE(numout,*) "                 in the tides, and making tidal analysis difficult." 
     139          ENDIF 
     140 
     141      ELSE         
     142          WRITE(numout,*) "       tides360: Gregorian calendar so using standard tides" 
     143      ENDIF 
     144 
    57145      CALL astronomic_angle 
    58146      CALL tide_pulse( pomega, ktide ,kc ) 
     
    69157      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2 
    70158      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac 
    71       !!---------------------------------------------------------------------- 
    72       ! 
    73       zqy = AINT( (nyear-1901.)/4. )        ! leap years since 1901 
    74       zsy = nyear - 1900.                   ! years since 1900 
    75       ! 
    76       zdj  = dayjul( nyear, nmonth, nday )  ! day number of year 
     159 
     160       
     161        ! JT 
     162        ! Tides are added as boundary conditions, and as tidal potential.  
     163        ! 
     164        ! For the Boundaries, the complex tide amplitudes are give for each point.  
     165        ! This gives the amplitude and the phase for each consititent.  
     166        ! The tidal frequency is calculated from Wave in tide.h90 via tide_pulse below. 
     167        ! A the start(ish) of everyday, astronomic_angle is called via tide_harmo  
     168        ! from SBC/sbctide.F90. 
     169        ! The astronomic_angle specifies the location of the moon and sun etc at the given  
     170        ! model time. these are used to update the tidal phase.  
     171        !  
     172        ! In the bdytides.F90 the function bdy_dta_tides calls 
     173        ! tide_init_elevation and tide_init_velocities (also in bdytides.F90)  
     174        ! every day. This uses the values from astro angles to update the phase of the  
     175        ! tidal constiuents as read in from the boundary files.  
     176        !  
     177        ! The tidal potential in (re)initialised every day in sbstide.F90 in the function 
     178        ! tide_init_potential. This uses the values from astro angles: 
     179        ! (v0tide + utide) and produces amp_pot and phi_pot. 
     180        ! These are then used in SBC/updtide.F90 (every timestep?) to set pot_astro. 
     181        ! 
     182        ! Both SBC/sbctide.F90 and bdy_dta_tides calculate zoff+z_arg which is the number of seconds since the beginning of the day. 
     183        ! the tidal phases are then corrected for this reset with the v0tide parameter, calucate by tide_vuf. 
     184        ! nodal correction is much smaller, with ftide (which affects the amplitude), and utide (which affects the phase). 
     185        !  
     186        !  
     187        ! As the phase of the tidal constituents for both the boundaries and the tidal potential 
     188        ! are adjusted by the astronomic_angle, we can adapt this one module to adapt the tides  
     189        ! for 360 day calendars. 
     190        !  
     191        ! There are different approaches to tides in a 360 day calendar. 
     192        !   1) (current), the tides are effectively reset to the first of the month. 
     193        !       therefore skip 31st's and repeat 29th and 30th of Feb 
     194        !       its happy with extra days of the month (doesn't crash for 30th Feb) 
     195        !       Tide is anchored to correct part of the year, but extra/missing days 
     196        !       are unrealistic, add noise to the system, and make least square tidal analysis fail 
     197        ! 
     198        !   2) Start the tides at the begining of the run and then let run continuously. 
     199        !       The tides drift throughout the year, so the equinox's are not at the correct part of the year. 
     200        !        
     201        !       This is the approach set up below 
     202        ! 
     203        !       2b) Adapt the equations to use decimal years (they sort of do, as they use day of year)  
     204        !           This would make the counting forward and backward from the origin easier  
     205        !           (the final step (going from DOY to mon and yr) would be removed) 
     206        ! 
     207        !   4) Adapt the equations that affect the location of the moon and tides.  
     208        !       This very likely to be a hugely complex job, that would affect the amphidromic systems, 
     209        !       As you're likely to need to change many/all of the tidal constants. this is then likely 
     210        !       to change the tidal frequencies, and so the the tidal wave speed, and hence the amphidromes,  
     211        !       co-tide and co-phase lines.  
     212        ! 
     213        !       This approach is not followed 
     214        ! 
     215        ! 
     216        ! To make the tide continueous for 360 and 365.25 day calendars, 
     217        ! firstly, I make temporary working yr/mn/day integers. 
     218        ! for the gregorian calendar these are simply set to nyear, nmonth and nday. 
     219        !  
     220        ! For a 360 day calendar, I count the days from 1900/1/1 to the current day  
     221        ! according to the the 360 day calendar. 
     222        ! I then count forward that many days according to the gregorian calendar. 
     223        ! therefore every 30 day month of the 360d model run, the tides move forward 30 days. 
     224 
     225 
     226 
     227    ! Questions: 
     228    ! Are the better ways of adding offets to dates in Nemo/Fortran? i.e. python timedelta from datetime 
     229    !       is there a leap year function in NEMO/Fortran? 
     230    ! Not sure if the algorithm is very stable for different origins. 
     231    !   nleap is corrected for 1900 not being a leap year. Needs updated for a different origin year 
     232    !   Does it work if its not starting on a leap year? 
     233    !   Does it work if its after the 28th Feb? 
     234    !   Does it work if the origin is after the start date of the run? 
     235    !   When adjusting the DOY and Y for the number of leap years, what happens if its more than 365 leap years? 
     236    !        
     237    ! h mean solar Longitude and s mean lunar Longitude and are functions of zhfrac, zday and zsy,  
     238    !   but the coeffiencents are not 1/86400:1:365 
     239    !   zday is the DOY corrected for the number of leap years since 1900.  
     240    !   So can run from 20 to 385. this wouldn't matter if the coefficients were 1/86400:1:365 
     241    !   should zsy and zday be updated so zday is between e.g. 1 and 365?? 
     242    ! 
     243    ! What are the impacts on the NWS if the tide drifts? 
     244    ! What is the impact on the NWS if the tide repeats/skips days? 
     245    !   can this make the model go unstable? 
     246     
     247    
     248 
     249    ! New variables defined for new code 
     250      INTEGER  ::   yr_org,mn_org,dy_org            !JT 
     251      REAL(wp) ::   sec_grg                         !JT 
     252      INTEGER  ::   yr_grg,mn_grg,dy_grg            !JT 
     253      INTEGER  ::   yr_360,mn_360,dy_360            !JT 
     254      INTEGER  ::   yr_wrk,mn_wrk,dy_wrk            !JT 
     255 
     256     
     257 
     258      LOGICAL  ::   ln_tide_drift_time_cont    ! Do we correct for a continueous time 
     259 
     260 
     261      ! INTEGER(KIND=8)  ::  days_since_origin  ! added to module 
     262      INTEGER  ::  init_yr, day_in_init_yr,nleap,init_doy 
     263      INTEGER  ::  init_doy_inc_l,yg_is_leap_mod,doy_grg 
     264      INTEGER,DIMENSION(12) ::  idayt, idays 
     265      INTEGER  ::   inc, ji 
     266      INTEGER  ::   ios 
     267 
     268 
     269      INTEGER  ::   yr_grg_2, mn_grg_2, dy_grg_2 
     270      REAL(wp) ::   sec_grg_2 
     271      REAL(wp) ::   fjulday_org       !: current julian day  
     272      ! REAL(wp) ::   days_since_origin_ymds2ju 
     273      INTEGER(KIND=8) ::   days_since_origin_ymds2ju_int 
     274 
     275      REAL(wp) ::   current_one_year 
     276      REAL(wp) ::   tmpju 
     277 
     278 
     279      REAL(wp) ::   jul_org_greg,jul_org_360,jul_pres_360 
     280          
     281      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 
     282       
     283       
     284      ! Currently hardcode the verbosity and the of the code 
     285      ! how to I read the calendar type 
     286      !ln_tide_360_cal = .TRUE. 
     287      !IF ( nleapy == 30)  ln_tide_360_cal = .TRUE. 
     288 
     289 
     290 
     291      !! Nameslist values 
     292 
     293 
     294      IF (ln_astro_verbose .AND. lwp) THEN 
     295        WRITE(numout,*) 'astro ' 
     296        WRITE(numout,*) 'astro -------------------------------------------------' 
     297      ENDIF 
     298       
     299 
     300      ln_tide_drift_time_cont = .FALSE. ! the same the original code 
     301 
     302      IF( nleapy == 30 ) THEN 
     303        IF ( ln_tide_drift ) THEN                
     304            ln_tide_drift_time_cont = .TRUE. 
     305        ENDIF 
     306        IF ( ln_tide_compress ) THEN                
     307            ! ################################################################## 
     308            ! ################################################################## 
     309            ! ################################################################## 
     310            ! ##### 
     311            ! #####  For the 360 day tide constituents,  
     312            ! #####  We only use days_since_origin for v0tide in tide_vuf. 
     313            ! #####   
     314            ! #####  To use the correct tide nodal correction (utide) 
     315            ! #####    (which is a small ajustment) 
     316            ! #####  use keep that linked to the gregorian dates. 
     317            ! ##### 
     318            ! ##### 
     319            ! ##### 
     320            ! #####  Therefore, we keep yr_wrk, mn_wrk, dy_wrk to equal 
     321            ! #####    nyear, nmonth, nday 
     322            ! ##### 
     323            ! ################################################################## 
     324            ! ################################################################## 
     325            ! ################################################################## 
     326 
     327            ln_tide_drift_time_cont = .FALSE. 
     328 
     329            ! ################################################################## 
     330            ! ################################################################## 
     331            ! ################################################################## 
     332            ! ##### 
     333            ! #####  NEMO4.0.4 
     334            ! #####  BUT!!! need to calc days_since_origin, so need to  
     335            ! #####  set ln_tide_drift_time_cont too true, then reset  
     336            ! #####  yr_wrk, mn_wrk, dy_wrk to equal nyear, nmonth, nday 
     337            ! ##### 
     338            ! ##### 
     339            ! ################################################################## 
     340            ! ################################################################## 
     341            ! ################################################################## 
     342            ln_tide_drift_time_cont = .TRUE. 
     343        ENDIF 
     344 
     345      ELSE 
     346        ln_tide_drift_time_cont = .FALSE. 
     347      ENDIF 
     348 
     349 
     350      IF (ln_astro_verbose .AND. lwp) THEN 
     351        WRITE(numout,*) 'astro ln_tide_drift_time_cont = ',ln_tide_drift_time_cont 
     352      ENDIF 
     353 
     354      !IF( ln_tide_360_cal ) THEN 
     355      !IF( nleapy == 30 ) THEN 
     356      IF( ln_tide_drift_time_cont ) THEN 
     357 
     358 
     359        ! clear and set current dates.  
     360        yr_360 = nyear 
     361        mn_360 = nmonth 
     362        dy_360 = nday 
     363 
     364 
     365        yr_grg = 0 
     366        mn_grg = 0 
     367        dy_grg = 0 
     368 
     369        yr_wrk = 0 
     370        mn_wrk = 0 
     371        dy_wrk = 0 
     372 
     373        ! Set the origin in the name list 
     374 
     375        yr_org = nn_tide_orig_yr 
     376        mn_org = nn_tide_orig_mn 
     377        dy_org = nn_tide_orig_dy 
     378 
     379         
     380        !IF (ln_tide_drift_time_cont_manual) THEN 
     381 
     382 
     383 
     384!            IF (ln_astro_verbose .AND. lwp) THEN 
     385!                WRITE(numout,*) 'astro: yr_360,yr_org,((yr_360-yr_org)*360)', yr_360,yr_org,((yr_360-yr_org)*360) 
     386!                WRITE(numout,*) 'astro: mn_360,mn_org,((mn_360-mn_org)*30)', mn_360,mn_org,((mn_360-mn_org)*30) 
     387!                WRITE(numout,*) 'astro: dy_360,dy_org,(dy_360-dy_org)', dy_360,dy_org,(dy_360-dy_org) 
     388!            ENDIF 
     389!             
     390!            ! how many days from 1900 in the 360 day calendar 
     391!            days_since_origin = ((yr_360-yr_org)*360) + ((mn_360-mn_org)*30) + (dy_360-dy_org) 
     392! 
     393!            ! first guess of what year this would be for the same numbers of days from 1/1/1900 in a gregorian calendar 
     394!            init_yr = yr_org + days_since_origin/365 
     395! 
     396!            ! was the initial estimated year a leap year? how many days in this year? 
     397!            day_in_init_yr = 365 
     398!            if (MOD(init_yr,4) == 0) day_in_init_yr = 366 
     399! 
     400! 
     401! 
     402!            !CALL ymds2ju_JT (yr_org, mn_org, dy_org, 0.0, fjulday_org,360.) 
     403! 
     404!            !IF (ln_astro_verbose) THEN 
     405!            !  IF(lwp) THEN 
     406!            !    WRITE(numout,*) 'astro: ymds2ju_JT yr_org, mn_org, dy_org,fjulday_org', yr_org, mn_org, dy_org,fjulday_org 
     407!            !  ENDIF 
     408!            !ENDIF 
     409! 
     410! 
     411!            !CALL ymds2ju( yr_org, mn_org, dy_org, 0.0, fjulday_org )  ! we assume that we start run at 00:00 
     412!            !IF( ABS(fjulday_org - REAL(NINT(fjulday_org),wp)) < 0.1 / rday )   fjulday_org = REAL(NINT(fjulday_org),wp)   ! avoid truncation error 
     413!            !fjulday_org = fjulday_org + 1.                             ! move back to the day at nit000 (and not at nit000 - 1) 
     414! 
     415!            !days_since_origin_ymds2ju_int = AINT(fjulday - fjulday_org) 
     416! 
     417!            IF (ln_astro_verbose .AND. lwp) THEN              
     418!                WRITE(numout,*) 'astro: days_since_origin,init_yr,day_in_init_yr', days_since_origin,init_yr,day_in_init_yr 
     419!                !WRITE(numout,*) 'astro: fjulday_org', fjulday_org 
     420!                !WRITE(numout,*) 'astro: fjulday', fjulday 
     421!                !WRITE(numout,*) 'astro: fjulday - fjulday_org', fjulday - fjulday_org 
     422!                !WRITE(numout,*) 'astro: days_since_origin_ymds2ju_int', days_since_origin_ymds2ju_int 
     423!            ENDIF 
     424! 
     425! 
     426!            ! how many leap years since the origin.  
     427!            nleap = (yr_360-yr_org)/4 - 1 !1900 is not a leap year 
     428!             
     429!            ! initial estimate of the day of year 
     430!            init_doy = MOD(days_since_origin,365) 
     431!             
     432!            ! correct the initial estimate for the DOY for the number of leap days since the origin 
     433!            init_doy_inc_l = init_doy - nleap 
     434! 
     435! 
     436!            IF (ln_astro_verbose .AND. lwp) THEN 
     437!                WRITE(numout,*) 'astro: nleap,init_doy,init_doy_inc_l',nleap,init_doy,init_doy_inc_l 
     438!            ENDIF 
     439!             
     440! 
     441!            ! The number of leap days could pull the  DOY before 0. 
     442!            ! in which case decrement the year, and reset the DOY. 
     443!            ! of the origin is 365 leap years ago, and initial DOY could be adjusted by more than one year.. 
     444!            ! Unlikely to be a prob, but need to remember if planning very long control runs. Need to think about this. 
     445! 
     446!            IF (init_doy_inc_l .LT. 0) THEN 
     447!                init_doy_inc_l = init_doy_inc_l+365 
     448!                init_yr = init_yr - 1  
     449!                IF (MOD(init_yr, 4) == 0 ) THEN 
     450!                    init_doy_inc_l = init_doy_inc_l + 1 
     451!                ENDIF 
     452!            ENDIF 
     453! 
     454!             
     455!            ! This gives the year and the day of year in the gregorian calendar 
     456!            yr_grg = init_yr     
     457!            doy_grg = init_doy_inc_l 
     458!            yg_is_leap_mod = MOD(yr_grg, 4) 
     459! 
     460!            IF (ln_astro_verbose .AND. lwp) THEN 
     461!                WRITE(numout,*) 'astro: yr_grg,doy_grg,yg_is_leap_mod',yr_grg,doy_grg,yg_is_leap_mod 
     462!            ENDIF 
     463! 
     464! 
     465!            ! Convert from day of year to month and day in the gregorian calendar. 
     466!            !   dayjul code adapted 
     467!            !   this perhaps should be a function, but not sure how to write one 
     468!            !   there may be this code functionality elsewhere in NEMO 
     469!            !!---------------------------------------------------------------------- 
     470!             
     471! 
     472!            ! what is the DOY of the first day of the month for each month. 
     473!            !   correct for leap years. 
     474!             
     475!            idays(1) = 0. 
     476!            idays(2) = 31. 
     477!            inc = 0. 
     478!            IF( yg_is_leap_mod == 0.)   inc = 1. 
     479! 
     480!            DO ji = 3, 12 
     481!                idays(ji)=idayt(ji)+inc 
     482!            END DO 
     483!         
     484!            ! cycle through the months. 
     485!            !   if the DOY is greater than the DOY of the first Day of Month 
     486!            !       Note the month. Calculate day of month by subtraction. 
     487!            !   Once beyond the correct month, the if statement won't be true, so wont calculate. 
     488! 
     489!            DO ji = 1, 12 
     490!                IF ( doy_grg .GE. idays(ji) )  THEN 
     491!                    mn_grg = ji 
     492!                    dy_grg = doy_grg-idays(ji) +1 
     493!                ENDIF 
     494!            END DO 
     495! 
     496! 
     497! 
     498! 
     499! 
     500!            IF(ln_astro_verbose .AND. lwp) THEN 
     501!                WRITE(numout,*) 'astro: mn_grg,dy_grg',mn_grg,dy_grg 
     502!                WRITE(numout,*) ' ' 
     503!                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_360,mn_360,dy_360,yr_grg,mn_grg,dy_grg,doy_grg =',yr_360,mn_360,dy_360,yr_grg,mn_grg,dy_grg,doy_grg 
     504! 
     505!                WRITE(numout,*) ' ' 
     506!            ENDIF 
     507!             
     508! 
     509!             
     510!            IF (ln_astro_verbose .AND. lwp)  WRITE(numout,*) 'tide_mod_astro_ang_meth_1,',yr_grg, mn_grg, dy_grg 
     511 
     512 
     513        !ELSE ! ln_tide_drift_time_cont_manual 
     514             
     515             
     516            ! number of days since 15th October 1582, for namelist origin, in both calendars, and for current model day. 
     517             
     518            CALL ymds2ju_JT( yr_org,mn_org,dy_org, 0. ,jul_org_greg,365.24 ) 
     519            CALL ymds2ju_JT( yr_org,mn_org,dy_org, 0. ,jul_org_360,360. ) 
     520            CALL ymds2ju_JT( yr_360,mn_360,dy_360, 0. ,jul_pres_360,360. ) 
     521 
     522            ! Calculate the days since the origin: days_since_origin_ymds2ju_int 
     523            ! How many days between the current day, and the origin, in the 360 day calendar. 
     524            days_since_origin_ymds2ju_int = jul_pres_360 - jul_org_360 
     525 
     526            IF (ln_astro_verbose .AND. lwp) THEN 
     527                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : jul_org_360,jul_pres_360,jul_pres_360 - jul_org_360 =',jul_org_360,jul_pres_360,jul_pres_360 - jul_org_360 
     528                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : days_since_origin_ymds2ju_int, days_since_origin_ymds2ju_int mod 360 =',days_since_origin_ymds2ju_int,MOD( days_since_origin_ymds2ju_int ,360 ) 
     529                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_org,mn_org,dy_org, jul_org_greg =',yr_org,mn_org,dy_org, jul_org_greg 
     530            ENDIF 
     531 
     532            !add days_since_origin_ymds2ju_int days to the origin in the gregorian calendar. 
     533            CALL ju2ymds_JT( days_since_origin_ymds2ju_int + jul_org_greg, yr_grg, mn_grg, dy_grg, sec_grg,365.24 ) 
     534 
     535            IF (ln_astro_verbose .AND. lwp) THEN 
     536                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_grg, mn_grg, dy_grg =',yr_grg, mn_grg, dy_grg 
     537                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_360, mn_360, dy_360 =',yr_360, mn_360, dy_360 
     538                WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_org, mn_org, dy_org =',yr_org, mn_org, dy_org 
     539            ENDIF 
     540 
     541 
     542 
     543             
     544            IF (ln_astro_verbose .AND. lwp)  WRITE(numout,*) 'tide_mod_astro_ang_meth_2,',yr_grg, mn_grg, dy_grg 
     545 
     546        !ENDIF !ln_tide_drift_time_cont_manual 
     547 
     548        ! for 360 calendars, work with the pseudo gregorian dates 
     549        yr_wrk = yr_grg 
     550        mn_wrk = mn_grg 
     551        dy_wrk = dy_grg 
     552 
     553        days_since_origin = days_since_origin_ymds2ju_int 
     554 
     555         
     556        IF (ln_tide_compress) THEN         
     557            yr_wrk = nyear 
     558            mn_wrk = nmonth 
     559            dy_wrk = nday 
     560        ENDIF 
     561 
     562      ELSE 
     563 
     564        ! for gregorian calendars, work with the model gregorian dates 
     565        yr_wrk = nyear 
     566        mn_wrk = nmonth 
     567        dy_wrk = nday 
     568 
     569      ENDIF 
     570 
     571      ! continue with original code, using working year, month and day. 
     572 
     573      ! 
     574      zqy = AINT( (yr_wrk-1901.)/4. )        ! leap years since 1901 
     575      zsy = yr_wrk - 1900.                   ! years since 1900 
     576      ! 
     577      zdj  = dayjul( yr_wrk, mn_wrk, dy_wrk )  ! day number of year 
    77578      zday = zdj + zqy - 1.                 ! day number of year + No of leap yrs 
    78579                                            ! i.e. what would doy if every year = 365 day?? 
    79580      ! 
    80581      zhfrac = nsec_day / 3600.             ! The seconds of the day/3600 
     582 
     583        
    81584      ! 
    82585      !---------------------------------------------------------------------- 
     
    105608      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 
    106609 
     610 
     611 
     612      IF(ln_astro_verbose .AND. lwp) THEN 
     613          WRITE(numout,*) 
     614          WRITE(numout,*) 'tide_mod_astro_ang : yr_wrk,mn_wrk,dy_wrk=',yr_wrk,mn_wrk,dy_wrk 
     615          WRITE(numout,*) 'tide_mod_astro_ang : nyear, nmonth, nday,nsec_day=',nyear, nmonth, nday,nsec_day 
     616          WRITE(numout,*) 'tide_mod_astro_ang : sh_N,sh_T,sh_h,sh_s,sh_p1,sh_p=', sh_N,sh_T,sh_h,sh_s,sh_p1,sh_p 
     617          WRITE(numout,*) 'tide_mod_astro_ang : zsy,zday,zhfrac,rad=', zsy,zday,zhfrac,rad 
     618          WRITE(numout,*) 'tide_mod_astro_ang : zqy ,zdj,yr_wrk, mn_wrk, dy_wrk =', zqy ,zdj,yr_wrk, mn_wrk, dy_wrk 
     619          WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     620      ENDIF 
     621 
     622 
     623 
    107624      sh_N = MOD( sh_N ,2*rpi ) 
    108625      sh_s = MOD( sh_s ,2*rpi ) 
     
    166683      INTEGER  ::   jh 
    167684      REAL(wp) ::   zscale 
    168       REAL(wp) ::   zomega_T =  13149000.0_wp 
    169       REAL(wp) ::   zomega_s =    481267.892_wp 
    170       REAL(wp) ::   zomega_h =     36000.76892_wp 
    171       REAL(wp) ::   zomega_p =      4069.0322056_wp 
    172       REAL(wp) ::   zomega_n =      1934.1423972_wp 
    173       REAL(wp) ::   zomega_p1=         1.719175_wp 
    174       !!---------------------------------------------------------------------- 
    175       ! 
    176       zscale =  rad / ( 36525._wp * 86400._wp )  
     685      REAL(wp) ::   zomega_T =  13149000.0_wp                ! Mean Solar Day  ! degrees/century 
     686      REAL(wp) ::   zomega_s =    481267.892_wp              ! Sidereal Month  ! degrees/century 
     687      REAL(wp) ::   zomega_h  !=     36000.76892_wp           ! Tropical Year   ! degrees/century 
     688      REAL(wp) ::   zomega_p =      4069.0322056_wp          ! Moons Perigee   ! degrees/century 
     689      REAL(wp) ::   zomega_n =      1934.1423972_wp          ! Regression of Lunar Nodes  ! degrees/century 
     690      REAL(wp) ::   zomega_p1=         1.719175_wp           ! Perihelion      ! degrees/century 
     691      !!---------------------------------------------------------------------- 
     692 
     693      zomega_h = 36000.76892_wp                  ! Tropical Year   ! degrees/century 
     694      IF (( nleapy == 30 ) .AND. ln_tide_compress)      zomega_h = 36525.0_wp   ! 360 day Tropical Year   ! degrees/century (360 deg/360 days= 1deg/day. cent = 36525 
     695 
     696      ! 
     697      zscale =  rad / ( 36525._wp * 86400._wp ) ! Convert to radians per second.  
    177698      ! 
    178699      DO jh = 1, kc 
     
    183704            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale 
    184705      END DO 
     706     
     707      IF (ln_astro_verbose .AND. lwp) THEN 
     708 
     709          WRITE(numout,*) 'astro tide_pulse nleapy:',nleapy 
     710          WRITE(numout,*) 'astro tide_pulse zomega_h:',zomega_h 
     711          WRITE(numout,*) 'astro tide_pulse if zomega_h = 36000.76892 for 365.24 day year' 
     712          WRITE(numout,*) 'astro tide_pulse if zomega_h = 36525.00000 for 360.00 day year' 
     713 
     714  
     715          DO jh = 1, kc 
     716              WRITE(numout,*) 'astro tide_pulse const, pomega, period(hr):',Wave(ktide(jh))%cname_tide, pomega(jh),2*rpi/(3600.0_wp*pomega(jh)) 
     717          END DO 
     718 
     719      ENDIF 
    185720      ! 
    186721   END SUBROUTINE tide_pulse 
     
    200735      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents 
    201736      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   ! 
     737 
    202738      ! 
    203739      INTEGER ::   jh   ! dummy loop index 
    204740      !!---------------------------------------------------------------------- 
    205       ! 
     741 
     742      !! JT for compress 
     743      REAL(wp)                             :: hours_since_origin 
     744      REAL(wp), DIMENSION(kc) ::   pomega           ! pulsation in radians/s 
     745      REAL(wp), DIMENSION(kc) ::   freq_per_day,   v0linearslope,v0linearintercept             ! pulsation in radians/s  !offset,cycle_reset,freq,per_hr 
     746      !! JT for compress 
     747 
     748 
     749 
     750      IF ( ln_tide_compress) THEN 
     751 
     752        CALL tide_pulse( pomega, ktide ,kc ) 
     753 
     754        DO jh = 1, kc 
     755          !per_hr(jh) = (2*rpi/pomega(jh))/3600. 
     756          !freq(jh) = (2*rpi/per_hr(jh)) 
     757          !freq_per_day(jh) = freq(jh)*24 
     758          freq_per_day(jh) = pomega(jh) * 86400.0_wp 
     759          !cycle_reset(jh) = mod(hours_since_origin*freq(jh),2.*rpi) 
     760          v0linearslope(jh) =   - mod (   (-freq_per_day(jh)), (2*rpi)  ) 
     761          IF(ln_astro_verbose .AND. lwp) WRITE(numout,*) 'astro tide_vuf 1:',jh,kc,ktide(jh),v0linearslope(jh),freq_per_day(jh), pomega(jh),(2*rpi/pomega(jh))/3600.! * 86400.0_wp,freq(jh)*24,per_hr(jh),freq(jh) 
     762        ENDDO 
     763 
     764 
     765        !offset(1) = 0.10789890_wp 
     766        !offset(2) = 1.10897897_wp 
     767        !offset(3) = 2.11005903_wp 
     768        !offset(4) = 0.00000000_wp 
     769        !offset(5) = 3.47632710_wp 
     770        !offset(6) = 0.16751976_wp 
     771        !offset(7) = -0.05503165_wp 
     772        !offset(8) = 0.94604842_wp 
     773        !offset(9) = 6.10534877_wp 
     774        !offset(10) = 0.21579780_wp 
     775        !offset(11) = 0.00000000_wp 
     776        !offset(12) = 0.00000000_wp 
     777        !offset(13) = 0.00000000_wp 
     778        !offset(14) = 0.00000000_wp 
     779        !offset(15) = 3.14159265_wp 
     780        !offset(16) = 0.21833313_wp 
     781        !offset(17) = 5.50043837_wp 
     782        !offset(18) = 2.24841149_wp 
     783        !offset(19) = 0.01800173_wp 
     784 
     785        !v0linearintercept(1) = 0.11044027_wp 
     786        !v0linearintercept(2) = 1.11152799_wp 
     787        !v0linearintercept(3) = 2.11261570_wp 
     788        !v0linearintercept(4) = 0.00000000_wp 
     789        !v0linearintercept(5) = 3.49727335_wp 
     790        !v0linearintercept(6) = 0.17784035_wp 
     791        !v0linearintercept(7) = 6.21578523_wp 
     792        !v0linearintercept(8) = 0.93368764_wp 
     793        !v0linearintercept(9) = 6.10534496_wp 
     794        !v0linearintercept(10) = 0.22088055_wp 
     795        !v0linearintercept(11) = 0.00000000_wp 
     796        !v0linearintercept(12) = 0.00000000_wp 
     797        !v0linearintercept(13) = 0.00000000_wp 
     798        !v0linearintercept(14) = 0.00000000_wp 
     799        !v0linearintercept(15) = 3.14159265_wp 
     800 
     801        !v0linearintercept(1) = v0linearintercept(1) - 0.000000_wp 
     802        !v0linearintercept(2) = v0linearintercept(2) - 0.000000_wp 
     803        !v0linearintercept(3) = v0linearintercept(3) - 0_wp 
     804        !v0linearintercept(4) = v0linearintercept(4) - 0.165795_wp 
     805        !v0linearintercept(5) = v0linearintercept(5) + 2.821252_wp 
     806        !v0linearintercept(6) = v0linearintercept(6) + 0.479504_wp 
     807        !v0linearintercept(7) = v0linearintercept(7) - 2.175621_wp 
     808        !v0linearintercept(8) = v0linearintercept(8) + 1.900267_wp 
     809        !v0linearintercept(9) = v0linearintercept(9) + 0.107633_wp 
     810        !v0linearintercept(10) = v0linearintercept(10) - 0.000000_wp 
     811        !v0linearintercept(11) = v0linearintercept(11) - 0.000000_wp 
     812        !v0linearintercept(12) = v0linearintercept(12) - 0.225730_wp 
     813        !v0linearintercept(13) = v0linearintercept(13) - 0.238641_wp 
     814        !v0linearintercept(14) = v0linearintercept(14) - 3.005851_wp 
     815        !v0linearintercept(15) = v0linearintercept(15) - 0.000000_wp 
     816 
     817        !v0linearintercept(1) =   0.11044026999999999_wp 
     818        !v0linearintercept(2) =   1.11152798999999990_wp 
     819        !v0linearintercept(3) =   2.11261570000000010_wp 
     820        !v0linearintercept(4) =  -0.16579500000000000_wp 
     821        !v0linearintercept(5) =   6.31852534999999980_wp 
     822        !v0linearintercept(6) =   0.65734435000000002_wp 
     823        !v0linearintercept(7) =   4.04016423000000020_wp 
     824        !v0linearintercept(8) =   2.83395464000000000_wp 
     825        !v0linearintercept(9) =   6.21297795999999990_wp 
     826        !v0linearintercept(10) =  0.22088055000000001_wp 
     827        !v0linearintercept(11) =  0.00000000000000000_wp 
     828        !v0linearintercept(12) = -0.22572999999999999_wp 
     829        !v0linearintercept(13) = -0.23864099999999999_wp 
     830        !v0linearintercept(14) = -3.00585099999999980_wp 
     831        !v0linearintercept(15) =  3.14159265000000020_wp 
     832 
     833        v0linearintercept( 1) =   0.2208805500_wp   -  (rpi* 68.0_wp/180.0_wp) !   M2  1 
     834        v0linearintercept( 2) =   3.1186126191_wp  !   N2  2 
     835        v0linearintercept( 3) =   0.9305155436_wp  !  2N2  3 
     836        v0linearintercept( 4) =   0.0194858941_wp  !   S2  4 
     837        v0linearintercept( 5) =  -2.5213114949_wp  !   K2  5 
     838        v0linearintercept( 6) =   6.5970532125_wp  !   K1  6 
     839        v0linearintercept( 7) =   1.1115279900_wp  !   O1  7 
     840        v0linearintercept( 8) =   0.1104402700_wp  !   Q1  8 
     841        !     v0linearintercept( 9) =   4.2269096542_wp  !   P1  9 
     842        !v0linearintercept( 9) =  -2.0351042402_wp  !   P1  9  compress3 
     843        !v0linearintercept( 9) =  -2.0351042402_wp  - 2.6179938779914944 !   P1  9  compress4 
     844 
     845        v0linearintercept( 9) =   rpi* 345.0_wp/180.0_wp -  (rpi* 140.0_wp/180.0_wp) !   P1  9  compress4 
     846 
     847        v0linearintercept(10) =   3.1415926500_wp  !   M4 10 
     848        v0linearintercept(11) =   0.0000000000_wp  !   Mf 11 
     849        v0linearintercept(12) =   0.0000000000_wp  !   Mm 12 
     850        v0linearintercept(13) =   0.0000000000_wp  ! Msqm 13 
     851        v0linearintercept(14) =   0.0000000000_wp  !  Mtm 14 
     852        v0linearintercept(15) =  -0.0230244122_wp  !   S1 15 
     853        v0linearintercept(16) =   4.2565208698_wp  !  MU2 16 
     854        v0linearintercept(17) =   6.5001767059_wp  !  NU2 17 
     855        v0linearintercept(18) =   0.0000000000_wp    -  (rpi* 113.0_wp/180.0_wp) !   L2 18 
     856        v0linearintercept(19) =   0.0092971808_wp  !   T2 19  + rpi/2. 
     857 
     858        !v0linearintercept(1) = v0linearintercept(1) - 0.034975_wp    ! M2 
     859        !v0linearintercept(2) = v0linearintercept(2) - 0.030244_wp    ! N2 
     860        !v0linearintercept(3) = v0linearintercept(3) - 0.036046_wp    ! 2N2 
     861        !v0linearintercept(4) = v0linearintercept(4) + 0.002092_wp    ! S2 
     862        !v0linearintercept(5) = v0linearintercept(5) - 0.273826_wp    ! K2 
     863        !v0linearintercept(6) = v0linearintercept(6) - 0.144677_wp    ! K1 
     864        !v0linearintercept(7) = v0linearintercept(7) + 0.031938_wp    ! O1 
     865        !v0linearintercept(8) = v0linearintercept(8) - 0.812030_wp    ! Q1 
     866        !v0linearintercept(9) = v0linearintercept(9) + 2.109118_wp    ! P1 
     867        !v0linearintercept(10) = v0linearintercept(10) + 0.070021_wp    ! M4 
     868        !v0linearintercept(11) = v0linearintercept(11) - 0.000000_wp    ! Mf 
     869        !v0linearintercept(12) = v0linearintercept(12) - 0.000000_wp    ! Mm 
     870        !v0linearintercept(13) = v0linearintercept(13) - 0.000000_wp    ! Msqm 
     871        !v0linearintercept(14) = v0linearintercept(14) - 0.000000_wp    ! Mtm 
     872        !v0linearintercept(15) = v0linearintercept(15) - 0.035676_wp    ! S1 
     873        !v0linearintercept(16) = v0linearintercept(16) + 0.007598_wp    ! MU2 
     874        !v0linearintercept(17) = v0linearintercept(17) - 0.043060_wp    ! NU2 
     875        !v0linearintercept(18) = v0linearintercept(18) + 0.023561_wp    ! L2 
     876        !v0linearintercept(19) = v0linearintercept(19) + 0.025624_wp    ! T2 
     877 
     878        v0linearintercept(1) = v0linearintercept(1) - (rpi*2.003909_wp/180.0_wp)    ! M2 
     879        v0linearintercept(2) = v0linearintercept(2) - (rpi*1.732874_wp/180.0_wp)    ! N2 
     880        v0linearintercept(3) = v0linearintercept(3) - (rpi*2.065265_wp/180.0_wp)    ! 2N2 
     881        v0linearintercept(4) = v0linearintercept(4) + (rpi*0.119842_wp/180.0_wp)    ! S2 
     882        v0linearintercept(5) = v0linearintercept(5) - (rpi*15.689068_wp/180.0_wp)    ! K2 
     883        v0linearintercept(6) = v0linearintercept(6) - (rpi*8.289390_wp/180.0_wp)    ! K1 
     884        v0linearintercept(7) = v0linearintercept(7) + (rpi*1.829931_wp/180.0_wp)    ! O1 
     885        v0linearintercept(8) = v0linearintercept(8) - (rpi*46.525902_wp/180.0_wp)    ! Q1 
     886        v0linearintercept(9) = v0linearintercept(9) + (rpi*120.843575_wp/180.0_wp)    ! P1 
     887        v0linearintercept(10) = v0linearintercept(10) + (rpi*4.011896_wp/180.0_wp)    ! M4 
     888        v0linearintercept(11) = v0linearintercept(11) - (rpi*0.000000_wp/180.0_wp)    ! Mf 
     889        v0linearintercept(12) = v0linearintercept(12) - (rpi*0.000000_wp/180.0_wp)    ! Mm 
     890        v0linearintercept(13) = v0linearintercept(13) - (rpi*0.000000_wp/180.0_wp)    ! Msqm 
     891        v0linearintercept(14) = v0linearintercept(14) - (rpi*0.000000_wp/180.0_wp)    ! Mtm 
     892        v0linearintercept(15) = v0linearintercept(15) - (rpi*2.044069_wp/180.0_wp)    ! S1 
     893        v0linearintercept(16) = v0linearintercept(16) + (rpi*0.435315_wp/180.0_wp)    ! MU2 
     894        v0linearintercept(17) = v0linearintercept(17) - (rpi*2.467160_wp/180.0_wp)    ! NU2 
     895        v0linearintercept(18) = v0linearintercept(18) + (rpi*1.349939_wp/180.0_wp)    ! L2 
     896        v0linearintercept(19) = v0linearintercept(19) + (rpi*1.468170_wp/180.0_wp)    ! T2 
     897 
     898 
     899        ! wave data. 
     900 
     901        !Wave( 1) = tide(  'M2'     , 0.242297 ,    2   ,  2 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     902        !Wave( 2) = tide(  'N2'     , 0.046313 ,    2   ,  2 , -3 ,  2 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     903        !Wave( 3) = tide( '2N2'     , 0.006184 ,    2   ,  2 , -4 ,  2 ,  2 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    78   ) 
     904        !Wave( 4) = tide(  'S2'     , 0.113572 ,    2   ,  2 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,     0   ) 
     905        !Wave( 5) = tide(  'K2'     , 0.030875 ,    2   ,  2 ,  0 ,  2 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   , -2   , 0 ,   235   ) 
     906        !!              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     907        !Wave( 6) = tide(  'K1'     , 0.142408 ,    1   ,  1 ,  0 ,  1 ,  0 ,  0  ,  -90  ,  0   ,  0   , -1   ,  0   , 0 ,   227   ) 
     908        !Wave( 7) = tide(  'O1'     , 0.101266 ,    1   ,  1 , -2 ,  1 ,  0 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   ) 
     909        !Wave( 8) = tide(  'Q1'     , 0.019387 ,    1   ,  1 , -3 ,  1 ,  1 ,  0  ,  +90  ,  2   , -1   ,  0   ,  0   , 0 ,    75   ) 
     910        !Wave( 9) = tide(  'P1'     , 0.047129 ,    1   ,  1 ,  0 , -1 ,  0 ,  0  ,  +90  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
     911        !!              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     912        !Wave(10) = tide(  'M4'     , 0.000000 ,    4   ,  4 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    ) 
     913        !!              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     914        !Wave(11) = tide(  'Mf'     , 0.042017 ,    0   ,  0 ,  2 ,  0 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     915        !Wave(12) = tide(  'Mm'     , 0.022191 ,    0   ,  0 ,  1 ,  0 , -1 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,   73    ) 
     916        !Wave(13) = tide(  'Msqm'   , 0.000667 ,    0   ,  0 ,  4 , -2 ,  0 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     917        !Wave(14) = tide(  'Mtm'    , 0.008049 ,    0   ,  0 ,  3 ,  0 , -1 ,  0  ,    0  , -2   ,  0   ,  0   ,  0   , 0 ,   74    ) 
     918        !!              !           !          !        !    !    !    !    !     !       !      !      !      !      !   !         ! 
     919        !Wave(15) = tide(  'S1'     , 0.000000 ,    1   ,  1 ,  0 ,  0 ,  0 ,  0  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    )    
     920        !Wave(16) = tide(  'MU2'    , 0.005841 ,    2   ,  2 , -4 ,  4 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
     921        !Wave(17) = tide(  'NU2'    , 0.009094 ,    2   ,  2 , -3 ,  4 , -1 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,   78    )  
     922        !Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    ) 
     923        !Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
     924 
     925        !name list 
     926        !  clname(1)='Q1' 
     927        !  clname(2)='O1' 
     928        !  clname(3)='P1' 
     929        !  clname(4)='S1' 
     930        !  clname(5)='K1' 
     931        !  clname(6)='2N2' 
     932        !  clname(7)='MU2' 
     933        !  clname(8)='N2' 
     934        !  clname(9)='NU2' 
     935        !  clname(10)='M2' 
     936        !  clname(11)='L2' 
     937        !  clname(12)='T2' 
     938        !  clname(13)='S2' 
     939        !  clname(14)='K2' 
     940        !  clname(15)='M4' 
     941 
     942        ! ktide 8,7,9,15 
     943 
     944        !ktide =  
     945        !8 
     946        !7 
     947        !9 
     948        !15 
     949        !6 
     950        !3 
     951        !16 
     952        !2 
     953        !17 
     954        !1 
     955        !18 
     956        !19 
     957        !4 
     958        !5 
     959        !10 
     960 
     961 
     962 
     963 
     964 
     965 
     966 
     967 
     968 
     969 
     970 
     971 
     972        !NEMO4 
     973 
     974!        clname(1)='Q1' 
     975!        clname(2)='O1' 
     976!        clname(3)='P1' 
     977!        clname(4)='S1' 
     978!        clname(5)='K1' 
     979!        clname(6)='2N2' 
     980!        clname(7)='MU2' 
     981!        clname(8)='N2' 
     982!        clname(9)='NU2' 
     983!        clname(10)='M2' 
     984!        clname(11)='L2' 
     985!        clname(12)='T2' 
     986!        clname(13)='S2' 
     987!        clname(14)='K2' 
     988!        clname(15)='M4' 
     989!        ktide = [10,9,11,12,8,23,21,15,22,14,18,19,16,17,28] 
     990 
     991 
     992        v0linearintercept( 1) =   0.1104402700_wp  !   Q1  8 
     993        v0linearintercept( 2) =   1.1115279900_wp  !   O1  7 
     994        v0linearintercept( 3) =   rpi* 345.0_wp/180.0_wp -  (rpi* 140.0_wp/180.0_wp) !   P1  9  compress4 
     995        v0linearintercept( 4) =  -0.0230244122_wp  !   S1 15 
     996        v0linearintercept( 5) =   6.5970532125_wp  !   K1  6 
     997        v0linearintercept( 6) =   0.9305155436_wp  !  2N2  3 
     998        v0linearintercept( 7) =   4.2565208698_wp  !  MU2 16 
     999        v0linearintercept( 8) =   3.1186126191_wp  !   N2  2 
     1000        v0linearintercept( 9) =   6.5001767059_wp  !  NU2 17 
     1001        v0linearintercept(10) =   0.2208805500_wp   -  (rpi* 68.0_wp/180.0_wp) !   M2  1 
     1002        v0linearintercept(11) =   0.0000000000_wp    -  (rpi* 113.0_wp/180.0_wp) !   L2 18 
     1003        v0linearintercept(12) =   0.0092971808_wp  !   T2 19  + rpi/2. 
     1004        v0linearintercept(13) =   0.0194858941_wp  !   S2  4 
     1005        v0linearintercept(14) =  -2.5213114949_wp  !   K2  5 
     1006        v0linearintercept(15) =   3.1415926500_wp  !   M4 10 
     1007 
     1008 
     1009 
     1010        v0linearintercept( 1) = v0linearintercept( 1) - (rpi*46.525902_wp/180.0_wp)   ! Q1 
     1011        v0linearintercept( 2) = v0linearintercept( 2) + (rpi*1.829931_wp/180.0_wp)    ! O1 
     1012        v0linearintercept( 3) = v0linearintercept( 3) + (rpi*120.843575_wp/180.0_wp)  ! P1 
     1013        v0linearintercept( 4) = v0linearintercept( 4) - (rpi*2.044069_wp/180.0_wp)    ! S1 
     1014        v0linearintercept( 5) = v0linearintercept( 5) - (rpi*8.289390_wp/180.0_wp)    ! K1 
     1015        v0linearintercept( 6) = v0linearintercept( 6) - (rpi*2.065265_wp/180.0_wp)    ! 2N2 
     1016        v0linearintercept( 7) = v0linearintercept( 7) + (rpi*0.435315_wp/180.0_wp)    ! MU2 
     1017        v0linearintercept( 8) = v0linearintercept( 8) - (rpi*1.732874_wp/180.0_wp)    ! N2 
     1018        v0linearintercept( 9) = v0linearintercept( 9) - (rpi*2.467160_wp/180.0_wp)    ! NU2 
     1019        v0linearintercept(10) = v0linearintercept(10) - (rpi*2.003909_wp/180.0_wp)    ! M2 
     1020        v0linearintercept(11) = v0linearintercept(11) + (rpi*1.349939_wp/180.0_wp)    ! L2 
     1021        v0linearintercept(12) = v0linearintercept(12) + (rpi*1.468170_wp/180.0_wp)    ! T2 
     1022        v0linearintercept(13) = v0linearintercept(13) + (rpi*0.119842_wp/180.0_wp)    ! S2 
     1023        v0linearintercept(14) = v0linearintercept(14) - (rpi*15.689068_wp/180.0_wp)   ! K2 
     1024        v0linearintercept(14) = v0linearintercept(15) + (rpi*4.011896_wp/180.0_wp)    ! M4 
     1025 
     1026 
     1027 
     1028 
     1029 
     1030 
     1031        DO jh = 1, kc 
     1032          IF(ln_astro_verbose .AND. lwp) WRITE(numout,*) 'astro tide_vuf 2:',jh,days_since_origin,v0linearslope(jh),v0linearintercept(ktide(jh))!,cycle_reset(jh)! ,offset(jh)  
     1033        ENDDO 
     1034      ENDIF 
     1035 
     1036 
     1037 
     1038      !!---------------------------------------------------------------------- 
     1039      !!---------------------------------------------------------------------- 
     1040      !! JT 
     1041      !!!  phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     1042      !!---------------------------------------------------------------------- 
     1043      !!---------------------------------------------------------------------- 
     1044 
    2061045      DO jh = 1, kc 
    2071046         !  Phase of the tidal potential relative to the Greenwhich  
    2081047         !  meridian (e.g. the position of the fictuous celestial body). Units are radian: 
    209          pvt(jh) = sh_T * Wave( ktide(jh) )%nT    & 
    210             &    + sh_s * Wave( ktide(jh) )%ns    & 
    211             &    + sh_h * Wave( ktide(jh) )%nh    & 
    212             &    + sh_p * Wave( ktide(jh) )%np    & 
    213             &    + sh_p1* Wave( ktide(jh) )%np1   & 
    214             &    +        Wave( ktide(jh) )%shift * rad 
     1048         !  Linear with time 
     1049 
     1050         IF ( ln_tide_compress ) THEN 
     1051 
     1052             pvt(jh) = mod(  ( (v0linearslope(jh)*days_since_origin) + v0linearintercept(  ktide(jh)  ) ),  2*rpi)-(2*rpi) 
     1053         ELSE 
     1054             pvt(jh) = sh_T * Wave( ktide(jh) )%nT    & 
     1055                &    + sh_s * Wave( ktide(jh) )%ns    & 
     1056                &    + sh_h * Wave( ktide(jh) )%nh    & 
     1057                &    + sh_p * Wave( ktide(jh) )%np    & 
     1058                &    + sh_p1* Wave( ktide(jh) )%np1   & 
     1059                &    +        Wave( ktide(jh) )%shift * rad 
     1060         ENDIF 
    2151061         ! 
    2161062         !  Phase correction u due to nodal motion. Units are radian: 
     1063         !  Cyclical with time. Much smaller terms than pvt.  
    2171064         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   & 
    2181065            &    + sh_nu     * Wave( ktide(jh) )%nnu0   & 
     
    2231070         !  Nodal correction factor: 
    2241071         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula ) 
     1072 
     1073 
    2251074      END DO 
     1075 
     1076 
     1077      IF(ln_astro_verbose .AND. lwp) THEN 
     1078          DO jh = 1, kc 
     1079              WRITE(numout,*) 'astro tide_vuf 3:',jh,pvt(jh), put(jh), pcor(jh)       
     1080          END DO 
     1081      ENDIF 
     1082 
     1083 
     1084 
    2261085      ! 
    2271086   END SUBROUTINE tide_vuf 
     
    4161275   END FUNCTION dayjul 
    4171276 
     1277 
     1278 
     1279 
     1280 
     1281 
     1282    SUBROUTINE ju2ymds_JT (julian,year,month,day,sec,one_year) 
     1283    !--------------------------------------------------------------------- 
     1284    !  IMPLICIT NONE 
     1285    !- 
     1286    !  REAL,INTENT(IN) :: julian 
     1287    !- 
     1288    !  INTEGER,INTENT(OUT) :: year,month,day 
     1289    !  REAL,INTENT(OUT)    :: sec 
     1290    !- 
     1291    !  INTEGER :: julian_day 
     1292    !  REAL    :: julian_sec 
     1293    !--------------------------------------------------------------------- 
     1294    !  julian_day = INT(julian) 
     1295    !  julian_sec = (julian-julian_day)*one_day 
     1296    !- 
     1297    !  CALL ju2ymds_internal(julian_day,julian_sec,year,month,day,sec) 
     1298    !--------------------- 
     1299    !END SUBROUTINE ju2ymds 
     1300    !- 
     1301    !=== 
     1302    !- 
     1303    !SUBROUTINE ju2ymds_internal (julian_day,julian_sec,year,month,day,sec) 
     1304    !--------------------------------------------------------------------- 
     1305    !- This subroutine computes from the julian day the year, 
     1306    !- month, day and seconds 
     1307    !- 
     1308    !- In 1968 in a letter to the editor of Communications of the ACM 
     1309    !- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel 
     1310    !- and Thomas C. Van Flandern presented such an algorithm. 
     1311    !- 
     1312    !- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm 
     1313    !- 
     1314    !- In the case of the Gregorian calendar we have chosen to use 
     1315    !- the Lilian day numbers. This is the day counter which starts 
     1316    !- on the 15th October 1582. This is the day at which Pope 
     1317    !- Gregory XIII introduced the Gregorian calendar. 
     1318    !- Compared to the true Julian calendar, which starts some 7980 
     1319    !- years ago, the Lilian days are smaler and are dealt with easily 
     1320    !- on 32 bit machines. With the true Julian days you can only the 
     1321    !- fraction of the day in the real part to a precision of a 1/4 of 
     1322    !- a day with 32 bits. 
     1323    !--------------------------------------------------------------------- 
     1324      IMPLICIT NONE 
     1325 
     1326 
     1327    !- 
     1328       
     1329      REAL,INTENT(IN)    :: julian,one_year 
     1330    !- 
     1331      INTEGER,INTENT(OUT) :: year,month,day 
     1332      REAL,INTENT(OUT)    :: sec 
     1333    !- 
     1334      INTEGER :: l,n,i,jd,j,d,m,y,ml 
     1335      INTEGER :: add_day 
     1336      REAL :: eps_day 
     1337 
     1338      REAL,PARAMETER :: one_day = 86400.0 
     1339    !--------------------------------------------------------------------- 
     1340 
     1341      INTEGER :: julian_day 
     1342      REAL    :: julian_sec 
     1343      INTEGER :: mon_len(12) 
     1344    !--------------------------------------------------------------------- 
     1345     
     1346    IF  ( (one_year > 365.0).AND.(one_year < 366.0) )  THEN 
     1347      mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/) 
     1348    ELSE IF ( ABS(one_year-365.0) <= EPSILON(one_year) ) THEN 
     1349      mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/) 
     1350    ELSE IF ( ABS(one_year-366.0) <= EPSILON(one_year) ) THEN 
     1351      mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/) 
     1352    ELSE IF ( ABS(one_year-360.0) <= EPSILON(one_year) ) THEN 
     1353      mon_len(:)=(/30,30,30,30,30,30,30,30,30,30,30,30/) 
     1354    ENDIF 
     1355 
     1356 
     1357 
     1358      julian_day = INT(julian) 
     1359      julian_sec = (julian-julian_day)*one_day 
     1360 
     1361      eps_day = SPACING(one_day) 
     1362    !  lock_one_year = .TRUE. 
     1363    !- 
     1364      jd = julian_day 
     1365      sec = julian_sec 
     1366      IF (sec > (one_day-eps_day)) THEN 
     1367        add_day = INT(sec/one_day) 
     1368        sec = sec-add_day*one_day 
     1369        jd = jd+add_day 
     1370      ENDIF 
     1371      IF (sec < -eps_day) THEN 
     1372        sec = sec+one_day 
     1373        jd = jd-1 
     1374      ENDIF 
     1375    !- 
     1376      IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN 
     1377    !-- Gregorian 
     1378        jd = jd+2299160 
     1379    !- 
     1380        l = jd+68569 
     1381        n = (4*l)/146097 
     1382        l = l-(146097*n+3)/4 
     1383        i = (4000*(l+1))/1461001 
     1384        l = l-(1461*i)/4+31 
     1385        j = (80*l)/2447 
     1386        d = l-(2447*j)/80 
     1387        l = j/11 
     1388        m = j+2-(12*l) 
     1389        y = 100*(n-49)+i+l 
     1390      ELSE IF (    (ABS(one_year-365.0) <= EPSILON(one_year)) & 
     1391     &         .OR.(ABS(one_year-366.0) <= EPSILON(one_year)) ) THEN 
     1392    !-- No leap or All leap 
     1393        !if ( ABS(one_year-365.0) <= EPSILON(one_year) ) mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/) 
     1394        !if ( ABS(one_year-366.0) <= EPSILON(one_year) ) mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/) 
     1395        y = jd/NINT(one_year) 
     1396        l = jd-y*NINT(one_year) 
     1397        m = 1 
     1398        ml = 0 
     1399        DO WHILE (ml+mon_len(m) <= l) 
     1400          ml = ml+mon_len(m) 
     1401          m = m+1 
     1402        ENDDO 
     1403        d = l-ml+1 
     1404      ELSE 
     1405    !-- others 
     1406        ml = NINT(one_year/12.) 
     1407        y = jd/NINT(one_year) 
     1408        l = jd-y*NINT(one_year) 
     1409        m = (l/ml)+1 
     1410        d = l-(m-1)*ml+1 
     1411      ENDIF 
     1412    !- 
     1413      day = d 
     1414      month = m 
     1415      year = y 
     1416    !------------------------------ 
     1417 
     1418 
     1419    END SUBROUTINE ju2ymds_JT 
     1420 
     1421 
     1422 
     1423 
     1424 
     1425 
     1426 
     1427 
     1428    !- 
     1429    SUBROUTINE ymds2ju_JT (year,month,day,sec,julian,one_year) 
     1430    !--------------------------------------------------------------------- 
     1431    !  IMPLICIT NONE 
     1432    !- 
     1433    !  INTEGER,INTENT(IN) :: year,month,day 
     1434    !  REAL,INTENT(IN)    :: sec 
     1435    !- 
     1436    !  REAL,INTENT(OUT) :: julian 
     1437    !- 
     1438    !  INTEGER :: julian_day 
     1439    !  REAL    :: julian_sec 
     1440    !--------------------------------------------------------------------- 
     1441    !  CALL ymds2ju_internal (year,month,day,sec,julian_day,julian_sec) 
     1442    !- 
     1443    !--------------------- 
     1444    !END SUBROUTINE ymds2ju 
     1445    !- 
     1446    !=== 
     1447    !- 
     1448    !SUBROUTINE ymds2ju_internal (year,month,day,sec,julian_day,julian_sec) 
     1449    !--------------------------------------------------------------------- 
     1450    !- Converts year, month, day and seconds into a julian day 
     1451    !- 
     1452    !- In 1968 in a letter to the editor of Communications of the ACM 
     1453    !- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel 
     1454    !- and Thomas C. Van Flandern presented such an algorithm. 
     1455    !- 
     1456    !- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm 
     1457    !- 
     1458    !- In the case of the Gregorian calendar we have chosen to use 
     1459    !- the Lilian day numbers. This is the day counter which starts 
     1460    !- on the 15th October 1582. 
     1461    !- This is the day at which Pope Gregory XIII introduced the 
     1462    !- Gregorian calendar. 
     1463    !- Compared to the true Julian calendar, which starts some 
     1464    !- 7980 years ago, the Lilian days are smaler and are dealt with 
     1465    !- easily on 32 bit machines. With the true Julian days you can only 
     1466    !- the fraction of the day in the real part to a precision of 
     1467    !- a 1/4 of a day with 32 bits. 
     1468    !--------------------------------------------------------------------- 
     1469      IMPLICIT NONE 
     1470    !- 
     1471      INTEGER,INTENT(IN) :: year,month,day 
     1472      REAL,INTENT(IN)    :: sec 
     1473      REAL,INTENT(IN)    :: one_year 
     1474    !- 
     1475    !  INTEGER,INTENT(OUT) :: julian_day 
     1476    !  REAL,INTENT(OUT)    :: julian_sec 
     1477      REAL,INTENT(OUT) :: julian 
     1478      INTEGER          :: julian_day 
     1479      REAL             :: julian_sec 
     1480    !- 
     1481      INTEGER :: jd,m,y,d,ml 
     1482      REAL,PARAMETER :: one_day = 86400.0 
     1483 
     1484 
     1485          INTEGER :: mon_len(12)   
     1486       
     1487        IF  ( (one_year > 365.0).AND.(one_year < 366.0) )  THEN 
     1488          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/) 
     1489        ELSE IF ( ABS(one_year-365.0) <= EPSILON(one_year) ) THEN 
     1490          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/) 
     1491        ELSE IF ( ABS(one_year-366.0) <= EPSILON(one_year) ) THEN 
     1492          mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/) 
     1493        ELSE IF ( ABS(one_year-360.0) <= EPSILON(one_year) ) THEN 
     1494          mon_len(:)=(/30,30,30,30,30,30,30,30,30,30,30,30/) 
     1495        ENDIF 
     1496 
     1497 
     1498    !--------------------------------------------------------------------- 
     1499     ! lock_one_year = .TRUE. 
     1500    !- 
     1501 
     1502 
     1503      m = month 
     1504      y = year 
     1505      d = day 
     1506 
     1507 
     1508        !--------------------------------------------------------------------- 
     1509 
     1510 
     1511    !- 
     1512    !- We deduce the calendar from the length of the year as it 
     1513    !- is faster than an INDEX on the calendar variable. 
     1514    !- 
     1515      IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN 
     1516    !-- "Gregorian" 
     1517        jd = (1461*(y+4800+INT((m-14)/12)))/4 & 
     1518     &      +(367*(m-2-12*(INT((m-14)/12))))/12 & 
     1519     &      -(3*((y+4900+INT((m-14)/12))/100))/4 & 
     1520     &      +d-32075 
     1521        jd = jd-2299160 
     1522      ELSE IF (    (ABS(one_year-365.0) <= EPSILON(one_year))  & 
     1523     &         .OR.(ABS(one_year-366.0) <= EPSILON(one_year)) ) THEN 
     1524    !-- "No leap" or "All leap" 
     1525        ml = SUM(mon_len(1:m-1)) 
     1526        jd = y*NINT(one_year)+ml+(d-1) 
     1527      ELSE 
     1528    !-- Calendar with regular month 
     1529        ml = NINT(one_year/12.) 
     1530        jd = y*NINT(one_year)+(m-1)*ml+(d-1) 
     1531      ENDIF 
     1532    !- 
     1533      julian_day = jd 
     1534      julian_sec = sec 
     1535 
     1536      julian = julian_day+julian_sec/one_day 
     1537 
     1538    !------------------------------ 
     1539    END SUBROUTINE ymds2ju_JT 
     1540 
     1541 
     1542 
     1543 
     1544 
     1545 
     1546 
     1547 
    4181548   !!====================================================================== 
    4191549END MODULE tide_mod 
Note: See TracChangeset for help on using the changeset viewer.