source: trunk/LMDZ.MARS/libf/phymars/iniorbit.F @ 3064

Last change on this file since 3064 was 3040, checked in by jbclement, 22 months ago

Mars PCM:
The variable 'timeperi' (defined in "planete_h.F90" and computed in "iniorbit.F") is renamed into 'lsperi' and thus slightly changed to be coherent to the solar longitude of perihelion in radian. It can now be used out of the box by other subroutines/programs like the PEM.
JBC

File size: 2.2 KB
Line 
1      SUBROUTINE iniorbit
2     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
3      use planete_h, only: aphelie, periheli, year_day, peri_day,
4     &                     obliquit, unitastr, e_elips, p_elips, lsperi
5      use comcstfi_h, only: pi
6      IMPLICIT NONE
7
8!=======================================================================
9! Initialisation of orbital parameters (stored in planete_h module)
10!=======================================================================
11
12!   Arguments:
13!   ----------
14
15      REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq
16
17!   Local:
18!   ------
19
20      REAL zxref,zanom,zz,zx0,zdx
21      INTEGER iter
22
23!-----------------------------------------------------------------------
24
25      pi=2.*asin(1.)
26
27      aphelie =paphelie
28      periheli=pperiheli
29      year_day=pyear_day
30      obliquit=pobliq
31      peri_day=pperi_day
32
33      PRINT*,'iniorbit: Perihelion in Mkm  ',periheli
34      PRINT*,'iniorbit: Aphelion  in Mkm  ',aphelie
35      PRINT*,'iniorbit: Obliquity in degrees  :',obliquit
36      unitastr=149.597927 ! 1 UA, in Mkm
37      e_elips=(aphelie-periheli)/(periheli+aphelie)
38      p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
39
40      print*,'iniorbit: e_elips',e_elips
41      print*,'iniorbit: p_elips',p_elips
42
43!-----------------------------------------------------------------------
44! compute polar angle and distance to the Sun:
45! -------------------------------------------------------
46
47!  compute mean anomaly zanom
48
49      zz=(year_day-pperi_day)/year_day
50      zanom=2.*pi*(zz-nint(zz))
51      zxref=abs(zanom)
52      PRINT*,'iniorbit: zanom  ',zanom
53
54!  solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
55!  using Newton method
56
57      zx0=zxref+e_elips*sin(zxref)
58      DO iter=1,100
59         zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
60         if(abs(zdx).le.(1.e-12)) exit
61         zx0=zx0+zdx
62      ENDDO
63
64      zx0=zx0+zdx
65      if(zanom.lt.0.) zx0=-zx0
66      PRINT*,'iniorbit: zx0   ',zx0
67
68      lsperi=-2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
69      lsperi = modulo(lsperi,2.*pi)
70      PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=',
71     &       lsperi*180./pi
72
73      END
Note: See TracBrowser for help on using the repository browser.