Changeset 15490 for NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate
- Timestamp:
- 2021-11-10T12:33:18+01:00 (3 years ago)
- 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 6 6 !! History : 1.0 ! 2007 (O. Le Galloudec) Original code 7 7 !!---------------------------------------------------------------------- 8 !JT USE oce ! ocean dynamics and tracers 9 USE lib_mpp ! distributed memory computing library 10 8 11 USE dom_oce ! ocean space and time domain 9 12 USE phycst ! physical constant 10 13 USE daymod ! calendar 14 USE in_out_manager ! I/O units 15 USE ioipsl , ONLY : ymds2ju ! for calendar 16 11 17 12 18 IMPLICIT NONE … … 34 40 REAL(wp) :: sh_I, sh_x1ra, sh_N ! 35 41 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 36 53 !!---------------------------------------------------------------------- 37 54 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 47 64 48 65 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 49 69 !!---------------------------------------------------------------------- 50 70 !!---------------------------------------------------------------------- … … 55 75 !!---------------------------------------------------------------------- 56 76 ! 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 ) 92 901 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 ) 96 902 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 57 145 CALL astronomic_angle 58 146 CALL tide_pulse( pomega, ktide ,kc ) … … 69 157 REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2 70 158 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 77 578 zday = zdj + zqy - 1. ! day number of year + No of leap yrs 78 579 ! i.e. what would doy if every year = 365 day?? 79 580 ! 80 581 zhfrac = nsec_day / 3600. ! The seconds of the day/3600 582 583 81 584 ! 82 585 !---------------------------------------------------------------------- … … 105 608 sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad 106 609 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 107 624 sh_N = MOD( sh_N ,2*rpi ) 108 625 sh_s = MOD( sh_s ,2*rpi ) … … 166 683 INTEGER :: jh 167 684 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. 177 698 ! 178 699 DO jh = 1, kc … … 183 704 & + zomega_p1* Wave( ktide(jh) )%np1 ) * zscale 184 705 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 185 720 ! 186 721 END SUBROUTINE tide_pulse … … 200 735 INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents 201 736 REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor ! 737 202 738 ! 203 739 INTEGER :: jh ! dummy loop index 204 740 !!---------------------------------------------------------------------- 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 206 1045 DO jh = 1, kc 207 1046 ! Phase of the tidal potential relative to the Greenwhich 208 1047 ! 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 215 1061 ! 216 1062 ! Phase correction u due to nodal motion. Units are radian: 1063 ! Cyclical with time. Much smaller terms than pvt. 217 1064 put(jh) = sh_xi * Wave( ktide(jh) )%nksi & 218 1065 & + sh_nu * Wave( ktide(jh) )%nnu0 & … … 223 1070 ! Nodal correction factor: 224 1071 pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula ) 1072 1073 225 1074 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 226 1085 ! 227 1086 END SUBROUTINE tide_vuf … … 416 1275 END FUNCTION dayjul 417 1276 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 418 1548 !!====================================================================== 419 1549 END MODULE tide_mod
Note: See TracChangeset
for help on using the changeset viewer.
