source: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 @ 3002

Last change on this file since 3002 was 2985, checked in by romain.vande, 3 years ago

Generic PEM :

Adapt PEM to run with the generic model.
(CPP_STD keyword to exclude some part of the code at the compilation)

RV

File size: 4.1 KB
Line 
1      MODULE recomp_orb_param_mod
2!=======================================================================
3!
4! Purpose: Recompute orbit parameters based on values by Laskar et al.,
5!
6! Author: RV
7!=======================================================================
8      IMPLICIT NONE
9
10      CONTAINS
11
12      SUBROUTINE recomp_orb_param(final_iter)
13
14      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp
15#ifndef CPP_STD     
16      USE comconst_mod, ONLY: pi
17      USE planete_h, ONLY: e_elips, obliquit, timeperi, periheli,aphelie,p_elips,peri_day,year_day
18#else
19      use planete_mod, only: e_elips, obliquit, timeperi,periastr,apoastr,p_elips,peri_day,year_day
20      USE comcstfi_mod, only: pi
21
22#endif
23
24      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
25                                end_lask_param_mod,last_ilask
26
27      IMPLICIT NONE
28
29!--------------------------------------------------------
30! Input Variables
31!--------------------------------------------------------
32
33      integer,intent(in) :: final_iter ! Number of year iteration of the PEM
34
35!--------------------------------------------------------
36! Output Variables
37!--------------------------------------------------------
38
39!--------------------------------------------------------
40! Local variables
41!--------------------------------------------------------
42
43      real :: Year                      ! Year of the simulation
44      real :: timeperi_ls               ! time of peri in ls
45      integer ilask                     ! Loop variable
46      real :: halfaxe                   ! Million of km
47      real :: unitastr
48
49      real xa,xb,xc,ya,yb
50
51
52      ! **********************************************************************
53      ! 0. Initializations
54      ! **********************************************************************
55#ifdef CPP_STD
56      real aphelie
57      real periheli
58
59      aphelie=apoastr
60      periheli=periastr
61#endif
62
63          Year=year_bp_ini+year_PEM+final_iter
64
65          timeperi_ls=timeperi*360/2/pi
66
67!          Year=-953200
68          print *, "recomp_orb_param, Old year=", year_bp_ini+year_PEM
69          print *, "recomp_orb_param, New year=", Year
70          print *, "recomp_orb_param, Old obl=", obliquit
71          print *, "recomp_orb_param, Old ex=", e_elips
72          print *, "recomp_orb_param, Old lsp=", timeperi_ls
73          print *, "recomp_orb_param, Old timeperi=", timeperi
74
75        xc=Year
76
77        do ilask=last_ilask,1,-1
78
79          xa=yearlask(ilask)
80          xb=yearlask(ilask+1)
81
82           if(yearlask(ilask) .GT.Year) then
83             if(var_obl) then
84              ya=oblask(ilask)
85              yb=oblask(ilask+1)
86              obliquit=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
87             endif
88             if(var_ex) then
89              ya=exlask(ilask)
90              yb=exlask(ilask+1)
91              e_elips=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
92             endif
93             if(var_lsp) then
94               if(lsplask(ilask)-lsplask(ilask+1) .gt.200) then
95                 if(lsplask(ilask).lt.lsplask(ilask+1)) then
96                   lsplask(ilask)=lsplask(ilask)+360
97                 else
98                   lsplask(ilask+1)=lsplask(ilask+1)+360
99                 endif
100               endif
101
102              ya=lsplask(ilask)
103              yb=lsplask(ilask+1)
104              timeperi_ls=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb)
105              timeperi_ls=modulo(timeperi_ls,360.)
106              halfaxe=227.94
107              timeperi=timeperi_ls*2*pi/360
108              periheli = halfaxe*(1-e_elips)
109              aphelie =  halfaxe*(1+e_elips)
110              unitastr=149.597927
111              p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
112#ifndef CPP_STD 
113              call call_dayperi(timeperi,e_elips,peri_day,year_day)
114#endif
115             endif
116              exit
117           endif
118        enddo
119
120       print *, "recomp_orb_param, Final year of the PEM run:", Year
121       print *, "recomp_orb_param, New obl=", obliquit
122       print *, "recomp_orb_param, New ex=", e_elips
123       print *, "recomp_orb_param, New timeperi=", timeperi
124
125      END SUBROUTINE recomp_orb_param
126
127!********************************************************************************   
128     
129      END MODULE recomp_orb_param_mod
Note: See TracBrowser for help on using the repository browser.