Index: trunk/LMDZ.MARS/libf/phymars/call_dayperi.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/call_dayperi.F	(revision 2228)
+++ trunk/LMDZ.MARS/libf/phymars/call_dayperi.F	(revision 2228)
@@ -0,0 +1,46 @@
+c ********************************************************
+      subroutine call_dayperi(Lsperi,e_elips,dayperi,year_day)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Objet:
+c   ------
+c   Computing the Martian date (number of sols since Ls=0 = spring
+c   equinox) at perihelion
+c
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   Lsperi      solar longitude (Ls) of perohelion (rad)
+c   e_elips       Excentricity
+c
+c   output
+c   ------
+c   dayperi       Martian date at perihelion (sol)
+c   real year_day      ! number of sols per Mars yar
+c-----------------------------------------------------------------------
+
+c arguments:
+c ----------
+       real Lsperi,e_elips,dayperi
+       real year_day      ! number of sols per Mars yar
+
+
+c      Local
+       real x1, x2,pi
+
+
+       pi=2.*asin(1.)
+
+       x1 = sqrt( (1-e_elips)/(1+e_elips) )
+       x2 = e_elips*sqrt(1-e_elips**2)
+
+       dayperi = 0.5*(year_day/pi)*
+     &   ( 2*atan(x1*tan(0.5*Lsperi))
+     &     -x2*sin(Lsperi)/(1+e_elips*cos(Lsperi)) )
+       if(dayperi.lt.0) dayperi=dayperi+year_day
+       return
+       end
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2226)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2228)
@@ -122,4 +122,8 @@
       character(len=44) :: txt
 
+c   New flag to compute paleo orbital configurations + few variables JN
+      REAL halfaxe, excentric, Lsperi
+      Logical paleomars
+
 c=======================================================================
 
@@ -151,6 +155,8 @@
       periheli = 206.66          ! minimum sun-mars distance (Mkm) ~206.66
       aphelie = 249.22           ! maximum sun-mars distance (Mkm) ~249.22
+      halfaxe = 227.94           ! demi-grand axe de l'ellipse
       peri_day =  485.           ! perihelion date (sols since N. Spring)
       obliquit = 25.2            ! Obliquity (deg) ~25.2         
+      excentric = 0.0934         ! Eccentricity (0.0934)         
  
 c     Planetary Boundary Layer and Turbulence parameters 
@@ -464,19 +470,44 @@
 c Orbital parameters
 c ------------------
-      print *,'Min. distance Sun-Mars (Mkm)?'
-      call getin("periheli",periheli)
-      write(*,*) " periheli = ",periheli
-
-      print *,'Max. distance Sun-Mars (Mkm)?'
-      call getin("aphelie",aphelie)
-      write(*,*) " aphelie = ",aphelie
-
-      print *,'Day of perihelion?'
-      call getin("periday",peri_day)
-      write(*,*) " periday = ",peri_day
-
-      print *,'Obliquity?'
-      call getin("obliquit",obliquit)
-      write(*,*) " obliquit = ",obliquit
+      paleomars=.false. ! Default: no water ice reservoir
+      call getin("paleomars",paleomars)
+      if (paleomars==.true.) then
+        write(*,*) "paleomars=", paleomars
+        write(*,*) "Orbital parameters from callphys.def"
+        write(*,*) "Enter eccentricity & Lsperi"
+        print *, 'Martian eccentricity (0<e<1) ?'
+        call getin('excentric ',excentric)
+        write(*,*)"excentric =",excentric 
+        print *, 'Solar longitude of perihelion (0<Ls<360) ?'
+        call getin('Lsperi',Lsperi )
+        write(*,*)"Lsperi=",Lsperi 
+        Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
+        periheli = halfaxe*(1-excentric)
+        aphelie = halfaxe*(1+excentric)
+        call call_dayperi(Lsperi,excentric,peri_day,year_day)
+        write(*,*) "Corresponding orbital params for GCM"
+        write(*,*) " periheli = ",periheli
+        write(*,*) " aphelie = ",aphelie
+        write(*,*) "date of perihelion (sol)",peri_day
+      else
+        write(*,*) "paleomars=", paleomars
+        write(*,*) "Default present-day orbital parameters"
+        write(*,*) "Unless specified otherwise"
+        print *,'Min. distance Sun-Mars (Mkm)?'
+        call getin("periheli",periheli)
+        write(*,*) " periheli = ",periheli
+
+        print *,'Max. distance Sun-Mars (Mkm)?'
+        call getin("aphelie",aphelie)
+        write(*,*) " aphelie = ",aphelie
+
+        print *,'Day of perihelion?'
+        call getin("periday",peri_day)
+        write(*,*) " periday = ",peri_day
+
+        print *,'Obliquity?'
+        call getin("obliquit",obliquit)
+        write(*,*) " obliquit = ",obliquit
+      end if
  
 c  latitude/longitude
Index: trunk/LMDZ.MARS/libf/phymars/planete_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/planete_h.F90	(revision 2226)
+++ trunk/LMDZ.MARS/libf/phymars/planete_h.F90	(revision 2228)
@@ -2,5 +2,5 @@
       IMPLICIT NONE
 
-      REAL aphelie ! Aphelion, in Mkm
+      REAL aphelie  ! Aphelion, in Mkm
       REAL periheli ! Perihelion, in Mkm
       REAL year_day ! number of days in the year