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