source: trunk/LMDZ.MARS/libf/phymars/planete_h.F90

Last change on this file was 3095, checked in by emillour, 14 months ago

Mars PCM:
Some code tidying:
Made pi in module comcstfi_h a parameter (and not a variable recomputed at
various points by various routines) and added module routine init_comcstfi_h
to cleanly initialize module variables.
Moved iniorbit.F to be a module routine of planete_h since it initializes
(some of ) the module variables it contains.
EM

File size: 2.4 KB
Line 
1MODULE planete_h
2IMPLICIT NONE
3
4REAL,SAVE :: aphelie   ! Aphelion, in Mkm
5REAL,SAVE :: periheli  ! Perihelion, in Mkm
6REAL,SAVE :: year_day  ! Number of days in the year
7!$OMP THREADPRIVATE(aphelie,periheli,year_day)
8REAL,SAVE :: peri_day  ! Date of perihelion, in days
9REAL,SAVE :: obliquit  ! Obliquity of the planet, in degrees
10REAL,SAVE :: lmixmin
11!$OMP THREADPRIVATE(peri_day,obliquit,lmixmin)
12REAL,SAVE :: emin_turb
13REAL,SAVE :: coefvis
14REAL,SAVE :: coefir
15!$OMP THREADPRIVATE(emin_turb,coefvis,coefir)
16REAL,SAVE :: lsperi    ! Solar longitude of the perihelion, angle in rad
17REAL,SAVE :: e_elips   ! Orbit eccentricity
18REAL,SAVE :: p_elips   ! Ellipse semi-latus rectum
19!$OMP THREADPRIVATE(lsperi,e_elips,p_elips)
20REAL,PARAMETER :: unitastr=149.597927 ! Astronomical unit AU, in Mkm
21
22! NB: Ideally all module variables should be "protected"...
23
24CONTAINS
25     
26  SUBROUTINE iniorbit(paphelie,pperiheli,pyear_day,pperi_day,pobliq)
27  use comcstfi_h, only: pi
28  IMPLICIT NONE
29  ! initialize orbital parameters
30 
31  REAL,INTENT(IN) :: paphelie,pperiheli,pyear_day,pperi_day,pobliq
32
33  CHARACTER(LEN=8),PARAMETER :: myname="iniorbit"
34  REAL :: zxref,zanom,zz,zx0,zdx
35  INTEGER :: iter
36
37  ! copy over input values
38  aphelie =paphelie
39  periheli=pperiheli
40  year_day=pyear_day
41  obliquit=pobliq
42  peri_day=pperi_day
43 
44  write(*,*)myname,': Perihelion in Mkm  ',periheli
45  write(*,*)myname,': Aphelion  in Mkm  ',aphelie
46  write(*,*)myname,': Obliquity in degrees  :',obliquit
47     
48  ! compute orbit geometrical parameters
49  e_elips=(aphelie-periheli)/(periheli+aphelie)
50  p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
51     
52  write(*,*)myname,': e_elips',e_elips
53  write(*,*)myname,': p_elips',p_elips
54     
55  !  compute mean anomaly zanom
56
57  zz=(year_day-pperi_day)/year_day
58  zanom=2.*pi*(zz-nint(zz))
59  zxref=abs(zanom)
60  write(*,*)myname,': zanom  ',zanom
61
62  ! solve equation  zx0 - e * sin (zx0) = zxref for eccentric anomaly zx0
63  ! using Newton method
64
65  zx0=zxref+e_elips*sin(zxref)
66  DO iter=1,100
67    zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0))
68    if(abs(zdx).le.(1.e-12)) exit
69    zx0=zx0+zdx
70  ENDDO
71
72  zx0=zx0+zdx
73  if(zanom.lt.0.) zx0=-zx0
74  write(*,*)myname,': zx0   ',zx0
75
76  lsperi=-2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
77  lsperi = modulo(lsperi,2.*pi)
78  write(*,*)myname,': Perihelion solar long. Ls (deg)=', &
79                          lsperi*180./pi
80
81  END SUBROUTINE iniorbit
82
83END MODULE planete_h
Note: See TracBrowser for help on using the repository browser.