[1992] | 1 | SUBROUTINE solarlong(pday, psollong, pdist_sol) |
---|
[1201] | 2 | |
---|
[1992] | 3 | USE ioipsl |
---|
[5112] | 4 | USE lmdz_print_control, ONLY: lunout |
---|
[5142] | 5 | USE lmdz_planete, ONLY: aphelie, periheli, year_day, peri_day, obliquit, timeperi, e_elips, p_elips, unitastr |
---|
[5144] | 6 | USE lmdz_yomcst |
---|
[1201] | 7 | |
---|
[1992] | 8 | IMPLICIT NONE |
---|
[1201] | 9 | |
---|
[1992] | 10 | ! ======================================================================= |
---|
[1201] | 11 | |
---|
[1992] | 12 | ! Objet: |
---|
| 13 | ! ------ |
---|
[1201] | 14 | |
---|
[1992] | 15 | ! Calcul de la distance soleil-planete et de la declinaison |
---|
| 16 | ! en fonction du jour de l'annee. |
---|
[1201] | 17 | |
---|
| 18 | |
---|
[1992] | 19 | ! Methode: |
---|
| 20 | ! -------- |
---|
[1201] | 21 | |
---|
[1992] | 22 | ! Calcul complet de l'elipse |
---|
[1201] | 23 | |
---|
[1992] | 24 | ! Interface: |
---|
| 25 | ! ---------- |
---|
[1201] | 26 | |
---|
[1992] | 27 | ! Uncommon comprenant les parametres orbitaux. |
---|
[1201] | 28 | |
---|
[1992] | 29 | ! Arguments: |
---|
| 30 | ! ---------- |
---|
[1201] | 31 | |
---|
[1992] | 32 | ! Input: |
---|
| 33 | ! ------ |
---|
| 34 | ! pday jour de l'annee (le jour 0 correspondant a l'equinoxe) |
---|
| 35 | ! lwrite clef logique pour sorties de controle |
---|
[1201] | 36 | |
---|
[1992] | 37 | ! Output: |
---|
| 38 | ! ------- |
---|
| 39 | ! pdist_sol distance entre le soleil et la planete |
---|
| 40 | ! ( en unite astronomique pour utiliser la constante |
---|
| 41 | ! solaire terrestre 1370 Wm-2 ) |
---|
| 42 | ! pdecli declinaison ( en radians ) |
---|
[1201] | 43 | |
---|
[1992] | 44 | ! ======================================================================= |
---|
| 45 | ! ----------------------------------------------------------------------- |
---|
[1201] | 46 | |
---|
[1992] | 47 | ! arguments: |
---|
| 48 | ! ---------- |
---|
[1201] | 49 | |
---|
[1992] | 50 | REAL pday, pdist_sol, pdecli, psollong |
---|
| 51 | LOGICAL lwrite |
---|
[1201] | 52 | |
---|
[1992] | 53 | ! Local: |
---|
| 54 | ! ------ |
---|
[1201] | 55 | |
---|
[1992] | 56 | REAL zanom, xref, zx0, zdx, zteta, zz, pi |
---|
| 57 | INTEGER iter |
---|
| 58 | REAL :: pyear_day, pperi_day |
---|
| 59 | REAL :: jd_eq, jd_peri |
---|
| 60 | LOGICAL, SAVE :: first = .TRUE. |
---|
| 61 | !$OMP THREADPRIVATE(first) |
---|
[1201] | 62 | |
---|
[1992] | 63 | ! ----------------------------------------------------------------------- |
---|
| 64 | ! calcul de l'angle polaire et de la distance au soleil : |
---|
| 65 | ! ------------------------------------------------------- |
---|
[1201] | 66 | |
---|
[1992] | 67 | ! Initialisation eventuelle: |
---|
| 68 | IF (first) THEN |
---|
| 69 | CALL ioget_calendar(pyear_day) |
---|
| 70 | CALL ymds2ju(2000, 3, 21, 0., jd_eq) |
---|
| 71 | CALL ymds2ju(2001, 1, 4, 0., jd_peri) |
---|
| 72 | pperi_day = jd_peri - jd_eq |
---|
| 73 | pperi_day = r_peri + 180. |
---|
| 74 | WRITE (lunout, *) ' Number of days in a year = ', pyear_day |
---|
[5103] | 75 | ! CALL iniorbit(249.22,206.66,669.,485.,25.2) |
---|
[1992] | 76 | CALL iniorbit(152.59, 146.61, pyear_day, pperi_day, r_incl) |
---|
| 77 | first = .FALSE. |
---|
| 78 | END IF |
---|
[1201] | 79 | |
---|
[1992] | 80 | ! calcul de l'zanomalie moyenne |
---|
[1201] | 81 | |
---|
[5142] | 82 | zz = (pday - peri_day) / year_day |
---|
| 83 | pi = 2. * asin(1.) |
---|
| 84 | zanom = 2. * pi * (zz - nint(zz)) |
---|
[1992] | 85 | xref = abs(zanom) |
---|
| 86 | |
---|
| 87 | ! resolution de l'equation horaire zx0 - e * sin (zx0) = xref |
---|
| 88 | ! methode de Newton |
---|
| 89 | |
---|
| 90 | ! zx0=xref+e_elips*sin(xref) |
---|
[5142] | 91 | zx0 = xref + r_ecc * sin(xref) |
---|
[1992] | 92 | DO iter = 1, 10 |
---|
| 93 | ! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) |
---|
[5142] | 94 | zdx = -(zx0 - r_ecc * sin(zx0) - xref) / (1. - r_ecc * cos(zx0)) |
---|
[1992] | 95 | IF (abs(zdx)<=(1.E-7)) GO TO 120 |
---|
| 96 | zx0 = zx0 + zdx |
---|
| 97 | END DO |
---|
[5142] | 98 | 120 CONTINUE |
---|
[1992] | 99 | zx0 = zx0 + zdx |
---|
| 100 | IF (zanom<0.) zx0 = -zx0 |
---|
| 101 | |
---|
| 102 | ! zteta est la longitude solaire |
---|
| 103 | |
---|
| 104 | ! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) |
---|
[5142] | 105 | zteta = 2. * atan(sqrt((1. + r_ecc) / (1. - r_ecc)) * tan(zx0 / 2.)) |
---|
[1992] | 106 | |
---|
| 107 | psollong = zteta - timeperi |
---|
| 108 | |
---|
[5142] | 109 | IF (psollong<0.) psollong = psollong + 2. * pi |
---|
| 110 | IF (psollong>2. * pi) psollong = psollong - 2. * pi |
---|
[1992] | 111 | |
---|
[5142] | 112 | psollong = psollong * 180. / pi |
---|
[1992] | 113 | |
---|
| 114 | ! distance soleil |
---|
| 115 | |
---|
[5142] | 116 | pdist_sol = (1 - r_ecc * r_ecc) / (1 + r_ecc * cos(pi / 180. * (psollong - & |
---|
| 117 | (r_peri + 180.0)))) |
---|
[1992] | 118 | ! pdist_sol = (1-e_elips*e_elips) |
---|
| 119 | ! & /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0)))) |
---|
| 120 | ! ----------------------------------------------------------------------- |
---|
| 121 | ! sorties eventuelles: |
---|
| 122 | ! --------------------- |
---|
| 123 | |
---|
| 124 | ! IF (lwrite) THEN |
---|
| 125 | ! PRINT*,'jour de l"annee :',pday |
---|
| 126 | ! PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol |
---|
| 127 | ! PRINT*,'declinaison (en degres) :',pdecli*180./pi |
---|
| 128 | ! ENDIF |
---|
| 129 | |
---|
| 130 | END SUBROUTINE solarlong |
---|