Changeset 12045
- Timestamp:
- 2019-12-03T18:24:37+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90
r12042 r12045 478 478 sh_p = MOD( sh_p, 2*rpi ) 479 479 480 cosI = 0.913694997 -0.035692561 *cos(sh_N)480 cosI = 0.913694997_wp - 0.035692561_wp * COS( sh_N ) 481 481 482 482 sh_I = ACOS( cosI ) … … 485 485 sh_tgn2 = tan(sh_N/2.0) 486 486 487 at1 =atan(1.01883*sh_tgn2)488 at2 =atan(0.64412*sh_tgn2)489 490 sh_xi =-at1-at2+sh_N491 492 IF( sh_N > rpi ) sh_xi=sh_xi-2.0*rpi487 at1 = ATAN( 1.01883_wp * sh_tgn2 ) 488 at2 = ATAN( 0.64412_wp * sh_tgn2 ) 489 490 sh_xi = sh_N - at1 - at2 491 492 IF( sh_N > rpi ) sh_xi = sh_xi - 2.0_wp * rpi 493 493 494 494 sh_nu = at1 - at2 495 495 496 !---------------------------------------------------------------------- 497 ! For constituents l2 k1 k2 498 !---------------------------------------------------------------------- 499 500 tgI2 = tan(sh_I/2.0) 501 P1 = sh_p-sh_xi 502 503 t2 = tgI2*tgI2 504 t4 = t2*t2 505 sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 ) 506 507 p = sin(2.0*P1) 508 q = 1.0/(6.0*t2)-cos(2.0*P1) 509 sh_R = atan(p/q) 510 511 p = sin(2.0*sh_I)*sin(sh_nu) 512 q = sin(2.0*sh_I)*cos(sh_nu)+0.3347 513 sh_nuprim = atan(p/q) 514 515 s2 = sin(sh_I)*sin(sh_I) 516 p = s2*sin(2.0*sh_nu) 517 q = s2*cos(2.0*sh_nu)+0.0727 518 sh_nusec = 0.5*atan(p/q) 496 ! For computation of tidal constituents L2 K1 K2 497 tgI2 = tan( sh_I / 2.0_wp ) 498 P1 = sh_p - sh_xi 499 ! 500 t2 = tgI2 * tgI2 501 t4 = t2 * t2 502 sh_x1ra = SQRT( 1.0 - 12.0 * t2 * COS( 2.0 * P1 ) + 36.0_wp * t4 ) 503 ! 504 p = SIN( 2.0_wp * P1 ) 505 q = 1.0_wp / ( 6.0_wp * t2 ) - COS( 2.0_wp * P1 ) 506 sh_R = ATAN( p / q ) 507 ! 508 p = SIN( 2.0_wp * sh_I ) * SIN( sh_nu ) 509 q = SIN( 2.0_wp * sh_I ) * COS( sh_nu ) + 0.3347_wp 510 sh_nuprim = ATAN( p / q ) 511 ! 512 s2 = SIN( sh_I ) * SIN( sh_I ) 513 p = s2 * SIN( 2.0_wp * sh_nu ) 514 q = s2 * COS( 2.0_wp * sh_nu ) + 0.0727_wp 515 sh_nusec = 0.5_wp * ATAN( p / q ) 519 516 ! 520 517 END SUBROUTINE astronomic_angle … … 595 592 RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf ) 596 593 !!---------------------------------------------------------------------- 597 !!---------------------------------------------------------------------- 598 INTEGER, INTENT(in) :: kformula 599 ! 600 REAL(wp) :: zf 601 REAL(wp) :: zs, zf1, zf2 602 CHARACTER(LEN=3) :: clformula 594 !! *** FUNCTION nodal_factort *** 595 !! 596 !! ** Purpose : Compute amplitude correction factors due to nodal motion 597 !! 598 !! ** Reference : 599 !! S58) Schureman, P. (1958): Manual of Harmonic Analysis and 600 !! Prediction of Tides (Revised (1940) Edition (Reprinted 1958 601 !! with corrections). Reprinted June 2001). U.S. Department of 602 !! Commerce, Coast and Geodetic Survey Special Publication 603 !! No. 98. Washington DC, United States Government Printing 604 !! Office. 317 pp. DOI: 10.25607/OBP-155. 605 !!---------------------------------------------------------------------- 606 INTEGER, INTENT(in) :: kformula 607 ! 608 REAL(wp) :: zf 609 REAL(wp) :: zs, zf1, zf2 610 CHARACTER(LEN=3) :: clformula 603 611 !!---------------------------------------------------------------------- 604 612 ! 605 613 SELECT CASE( kformula ) 606 614 ! 607 CASE( 0 ) ! == formule0, solar waves615 CASE( 0 ) ! Formula 0, solar waves 608 616 zf = 1.0 609 617 ! 610 CASE( 1 ) ! == formule1, compound waves (78 x 78)611 zf=nodal_factort( 78)618 CASE( 1 ) ! Formula 1, compound waves (78 x 78) 619 zf=nodal_factort( 78 ) 612 620 zf = zf * zf 613 621 ! 614 CASE ( 4 ) ! == formule4, compound waves (78 x 235)615 zf1 = nodal_factort( 78 )622 CASE ( 4 ) ! Formula 4, compound waves (78 x 235) 623 zf1 = nodal_factort( 78 ) 616 624 zf = nodal_factort(235) 617 625 zf = zf1 * zf 618 626 ! 619 CASE( 18 ) ! == formule18, compound waves (78 x 78 x 78 )620 zf1 = nodal_factort( 78)627 CASE( 18 ) ! Formula 18, compound waves (78 x 78 x 78 ) 628 zf1 = nodal_factort( 78 ) 621 629 zf = zf1 * zf1 * zf1 622 630 ! 623 CASE( 20 ) ! == formula 20, compound waves ( 78 x 78 x 78 x 78 )624 zf1 = nodal_factort( 78)631 CASE( 20 ) ! Formula 20, compound waves ( 78 x 78 x 78 x 78 ) 632 zf1 = nodal_factort( 78 ) 625 633 zf = zf1 * zf1 * zf1 * zf1 626 634 ! 627 CASE( 73 ) ! == formule 73628 zs = sin(sh_I)629 zf = ( 2./3.-zs*zs)/0.5021630 ! 631 CASE( 74 ) ! == formule 74632 zs = sin(sh_I)633 zf = zs * zs / 0.1578 634 ! 635 CASE( 75 ) ! == formule 75636 zs = cos(sh_I/2)637 zf = sin(sh_I) * zs * zs / 0.3800638 ! 639 CASE( 76 ) ! == formule 76640 zf = sin(2*sh_I) / 0.7214641 ! 642 CASE( 78 ) ! == formule 78643 zs = cos(sh_I/2)644 zf = zs * zs * zs * zs / 0.9154 645 ! 646 CASE( 149 ) ! == formule 149647 zs = cos(sh_I/2)648 zf = zs *zs*zs*zs*zs*zs / 0.8758649 ! 650 CASE( 215 ) ! == formule 215651 zs = cos(sh_I/2)652 zf = zs *zs*zs*zs / 0.9154* sh_x1ra653 ! 654 CASE( 227 ) ! == formule 227655 zs = sin(2*sh_I)656 zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006)657 ! 658 CASE ( 235 ) ! == formule 235659 zs = sin(sh_I)660 zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981)635 CASE( 73 ) ! Formula 73 of S58 636 zs = SIN( sh_I ) 637 zf = ( 2.0_wp / 3.0_wp - zs * zs ) / 0.5021_wp 638 ! 639 CASE( 74 ) ! Formula 74 of S58 640 zs = SIN(sh_I) 641 zf = zs * zs / 0.1578_wp 642 ! 643 CASE( 75 ) ! Formula 75 of S58 644 zs = COS( sh_I / 2.0_wp ) 645 zf = SIN( sh_I ) * zs * zs / 0.3800_wp 646 ! 647 CASE( 76 ) ! Formula 76 of S58 648 zf = SIN( 2.0_wp * sh_I ) / 0.7214_wp 649 ! 650 CASE( 78 ) ! Formula 78 of S58 651 zs = COS( sh_I/2 ) 652 zf = zs * zs * zs * zs / 0.9154_wp 653 ! 654 CASE( 149 ) ! Formula 149 of S58 655 zs = COS( sh_I/2 ) 656 zf = zs * zs * zs * zs * zs * zs / 0.8758_wp 657 ! 658 CASE( 215 ) ! Formula 215 of S58 with typo correction (0.9154 instead of 0.9145) 659 zs = COS( sh_I/2 ) 660 zf = zs * zs * zs * zs / 0.9154_wp * sh_x1ra 661 ! 662 CASE( 227 ) ! Formula 227 of S58 663 zs = SIN( 2.0_wp * sh_I ) 664 zf = SQRT( 0.8965_wp * zs * zs + 0.6001_wp * zs * COS( sh_nu ) + 0.1006_wp ) 665 ! 666 CASE ( 235 ) ! Formula 235 of S58 667 zs = SIN( sh_I ) 668 zf = SQRT( 19.0444_wp * zs * zs * zs * zs + 2.7702_wp * zs * zs * cos( 2.0_wp * sh_nu ) + 0.0981_wp ) 661 669 ! 662 670 CASE DEFAULT … … 670 678 FUNCTION dayjul( kyr, kmonth, kday ) 671 679 !!---------------------------------------------------------------------- 672 !! *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE) 673 !!---------------------------------------------------------------------- 674 INTEGER,INTENT(in) :: kyr, kmonth, kday 675 ! 676 INTEGER,DIMENSION(12) :: idayt, idays 677 INTEGER :: inc, ji 678 REAL(wp) :: dayjul, zyq 679 ! 680 DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./ 681 !!---------------------------------------------------------------------- 682 ! 683 idays(1) = 0. 684 idays(2) = 31. 685 inc = 0. 686 zyq = MOD( kyr-1900. , 4. ) 687 IF( zyq == 0.) inc = 1. 680 !! *** FUNCTION dayjul *** 681 !! 682 !! Purpose : compute the Julian day 683 !!---------------------------------------------------------------------- 684 INTEGER,INTENT(in) :: kyr, kmonth, kday 685 ! 686 INTEGER,DIMENSION(12) :: idayt = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /) 687 INTEGER,DIMENSION(12) :: idays 688 INTEGER :: inc, ji, zyq 689 REAL(wp) :: dayjul 690 !!---------------------------------------------------------------------- 691 ! 692 idays(1) = 0 693 idays(2) = 31 694 inc = 0.0_wp 695 zyq = MOD( kyr - 1900 , 4 ) 696 IF( zyq == 0 ) inc = 1 688 697 DO ji = 3, 12 689 idays(ji) =idayt(ji)+inc698 idays(ji) = idayt(ji) + inc 690 699 END DO 691 dayjul = idays(kmonth) + kday700 dayjul = REAL( idays(kmonth) + kday, KIND=wp ) 692 701 ! 693 702 END FUNCTION dayjul
Note: See TracChangeset
for help on using the changeset viewer.