[1992] | 1 | SUBROUTINE solarlong(pday, psollong, pdist_sol) |
---|
[1201] | 2 | |
---|
[1992] | 3 | USE ioipsl |
---|
[2311] | 4 | USE print_control_mod, ONLY: lunout |
---|
[1201] | 5 | |
---|
[5274] | 6 | USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO & |
---|
| 7 | , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA & |
---|
| 8 | , R_ecc, R_peri, R_incl & |
---|
| 9 | , RA, RG, R1SA & |
---|
| 10 | , RSIGMA & |
---|
| 11 | , R, RMD, RMV, RD, RV, RCPD & |
---|
| 12 | , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12 & |
---|
| 13 | , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w & |
---|
| 14 | , RCW, RCS & |
---|
| 15 | , RLVTT, RLSTT, RLMLT, RTT, RATM & |
---|
| 16 | , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS & |
---|
| 17 | , RALPD, RBETD, RGAMD |
---|
| 18 | IMPLICIT NONE |
---|
[1201] | 19 | |
---|
[1992] | 20 | ! ======================================================================= |
---|
[1201] | 21 | |
---|
[1992] | 22 | ! Objet: |
---|
| 23 | ! ------ |
---|
[1201] | 24 | |
---|
[1992] | 25 | ! Calcul de la distance soleil-planete et de la declinaison |
---|
| 26 | ! en fonction du jour de l'annee. |
---|
[1201] | 27 | |
---|
| 28 | |
---|
[1992] | 29 | ! Methode: |
---|
| 30 | ! -------- |
---|
[1201] | 31 | |
---|
[1992] | 32 | ! Calcul complet de l'elipse |
---|
[1201] | 33 | |
---|
[1992] | 34 | ! Interface: |
---|
| 35 | ! ---------- |
---|
[1201] | 36 | |
---|
[1992] | 37 | ! Uncommon comprenant les parametres orbitaux. |
---|
[1201] | 38 | |
---|
[1992] | 39 | ! Arguments: |
---|
| 40 | ! ---------- |
---|
[1201] | 41 | |
---|
[1992] | 42 | ! Input: |
---|
| 43 | ! ------ |
---|
| 44 | ! pday jour de l'annee (le jour 0 correspondant a l'equinoxe) |
---|
| 45 | ! lwrite clef logique pour sorties de controle |
---|
[1201] | 46 | |
---|
[1992] | 47 | ! Output: |
---|
| 48 | ! ------- |
---|
| 49 | ! pdist_sol distance entre le soleil et la planete |
---|
| 50 | ! ( en unite astronomique pour utiliser la constante |
---|
| 51 | ! solaire terrestre 1370 Wm-2 ) |
---|
| 52 | ! pdecli declinaison ( en radians ) |
---|
[1201] | 53 | |
---|
[1992] | 54 | ! ======================================================================= |
---|
| 55 | ! ----------------------------------------------------------------------- |
---|
| 56 | ! Declarations: |
---|
| 57 | ! ------------- |
---|
[1201] | 58 | |
---|
[1992] | 59 | include "planete.h" |
---|
[1201] | 60 | |
---|
[5274] | 61 | |
---|
[1992] | 62 | ! arguments: |
---|
| 63 | ! ---------- |
---|
[1201] | 64 | |
---|
[1992] | 65 | REAL pday, pdist_sol, pdecli, psollong |
---|
| 66 | LOGICAL lwrite |
---|
[1201] | 67 | |
---|
[1992] | 68 | ! Local: |
---|
| 69 | ! ------ |
---|
[1201] | 70 | |
---|
[1992] | 71 | REAL zanom, xref, zx0, zdx, zteta, zz, pi |
---|
| 72 | INTEGER iter |
---|
| 73 | REAL :: pyear_day, pperi_day |
---|
| 74 | REAL :: jd_eq, jd_peri |
---|
| 75 | LOGICAL, SAVE :: first = .TRUE. |
---|
| 76 | !$OMP THREADPRIVATE(first) |
---|
[1201] | 77 | |
---|
[1992] | 78 | ! ----------------------------------------------------------------------- |
---|
| 79 | ! calcul de l'angle polaire et de la distance au soleil : |
---|
| 80 | ! ------------------------------------------------------- |
---|
[1201] | 81 | |
---|
[1992] | 82 | ! Initialisation eventuelle: |
---|
| 83 | IF (first) THEN |
---|
| 84 | CALL ioget_calendar(pyear_day) |
---|
| 85 | CALL ymds2ju(2000, 3, 21, 0., jd_eq) |
---|
| 86 | CALL ymds2ju(2001, 1, 4, 0., jd_peri) |
---|
| 87 | pperi_day = jd_peri - jd_eq |
---|
| 88 | pperi_day = r_peri + 180. |
---|
| 89 | WRITE (lunout, *) ' Number of days in a year = ', pyear_day |
---|
| 90 | ! call iniorbit(249.22,206.66,669.,485.,25.2) |
---|
| 91 | CALL iniorbit(152.59, 146.61, pyear_day, pperi_day, r_incl) |
---|
| 92 | first = .FALSE. |
---|
| 93 | END IF |
---|
[1201] | 94 | |
---|
[1992] | 95 | ! calcul de l'zanomalie moyenne |
---|
[1201] | 96 | |
---|
[1992] | 97 | zz = (pday-peri_day)/year_day |
---|
| 98 | pi = 2.*asin(1.) |
---|
| 99 | zanom = 2.*pi*(zz-nint(zz)) |
---|
| 100 | xref = abs(zanom) |
---|
| 101 | |
---|
| 102 | ! resolution de l'equation horaire zx0 - e * sin (zx0) = xref |
---|
| 103 | ! methode de Newton |
---|
| 104 | |
---|
| 105 | ! zx0=xref+e_elips*sin(xref) |
---|
| 106 | zx0 = xref + r_ecc*sin(xref) |
---|
| 107 | DO iter = 1, 10 |
---|
| 108 | ! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) |
---|
| 109 | zdx = -(zx0-r_ecc*sin(zx0)-xref)/(1.-r_ecc*cos(zx0)) |
---|
| 110 | IF (abs(zdx)<=(1.E-7)) GO TO 120 |
---|
| 111 | zx0 = zx0 + zdx |
---|
| 112 | END DO |
---|
| 113 | 120 CONTINUE |
---|
| 114 | zx0 = zx0 + zdx |
---|
| 115 | IF (zanom<0.) zx0 = -zx0 |
---|
| 116 | |
---|
| 117 | ! zteta est la longitude solaire |
---|
| 118 | |
---|
| 119 | ! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) |
---|
| 120 | zteta = 2.*atan(sqrt((1.+r_ecc)/(1.-r_ecc))*tan(zx0/2.)) |
---|
| 121 | |
---|
| 122 | psollong = zteta - timeperi |
---|
| 123 | |
---|
| 124 | IF (psollong<0.) psollong = psollong + 2.*pi |
---|
| 125 | IF (psollong>2.*pi) psollong = psollong - 2.*pi |
---|
| 126 | |
---|
| 127 | psollong = psollong*180./pi |
---|
| 128 | |
---|
| 129 | ! distance soleil |
---|
| 130 | |
---|
| 131 | pdist_sol = (1-r_ecc*r_ecc)/(1+r_ecc*cos(pi/180.*(psollong- & |
---|
| 132 | (r_peri+180.0)))) |
---|
| 133 | ! pdist_sol = (1-e_elips*e_elips) |
---|
| 134 | ! & /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0)))) |
---|
| 135 | ! ----------------------------------------------------------------------- |
---|
| 136 | ! sorties eventuelles: |
---|
| 137 | ! --------------------- |
---|
| 138 | |
---|
| 139 | ! IF (lwrite) THEN |
---|
| 140 | ! PRINT*,'jour de l"annee :',pday |
---|
| 141 | ! PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol |
---|
| 142 | ! PRINT*,'declinaison (en degres) :',pdecli*180./pi |
---|
| 143 | ! ENDIF |
---|
| 144 | |
---|
| 145 | RETURN |
---|
| 146 | END SUBROUTINE solarlong |
---|