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

Last change on this file since 2266 was 1382, checked in by emillour, 10 years ago

Mars GCM:
Some cleanup in iniorbit.F
FF+EM

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