Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/hgardfou.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/hgardfou.F	(revision 1198)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/hgardfou.F	(revision 1201)
@@ -1,3 +1,4 @@
 !
+! $Id$
       SUBROUTINE hgardfou (t,tsol,text)
       use dimphy
@@ -12,4 +13,6 @@
       REAL t(klon,klev), tsol(klon,nbsrf)
       CHARACTER*(*) text
+      character (len=20) :: modname = 'hgardfou'
+      character (len=80) :: abort_message
 C
       INTEGER i, k, nsrf
@@ -124,6 +127,6 @@
 c
       IF (.NOT. ok) THEN
-         PRINT*, 'hgardfou s arrete ', text
-         CALL abort
+         abort_message= 'hgardfou s arrete '//text
+         CALL abort_gcm (modname,abort_message,1)
       ENDIF
 
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/iniorbit.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/iniorbit.F	(revision 1201)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/iniorbit.F	(revision 1201)
@@ -0,0 +1,104 @@
+      SUBROUTINE iniorbit
+     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteur:
+c   -------
+c     Frederic Hourdin      22 Fevrier 1991
+c
+c   Objet:
+c   ------
+c    Initialisation du sous programme orbite qui calcule
+c    a une date donnee de l'annee de duree year_day commencant
+c    a l'equinoxe de printemps et dont le perihelie se situe
+c    a la date peri_day, la distance au soleil et la declinaison.
+c
+c   Interface:
+c   ----------
+c   - Doit etre appele avant d'utiliser orbite.
+c   - initialise une partie du common planete.h
+c
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   aphelie       \   aphelie et perihelie de l'orbite
+c   periheli      /   en millions de kilometres.
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+#include "planete.h"
+#include "YOMCST.h"
+
+c   Arguments:
+c   ----------
+
+      REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq
+
+c   Local:
+c   ------
+
+      REAL zxref,zanom,zz,zx0,zdx, pi
+      INTEGER iter
+
+c-----------------------------------------------------------------------
+
+      pi=2.*asin(1.)
+
+      aphelie =paphelie
+      periheli=pperiheli
+      year_day=pyear_day
+      obliquit=pobliq
+      peri_day=pperi_day
+
+      PRINT*,'Perihelie en Mkm  ',periheli
+      PRINT*,'Aphelie  en Mkm   ',aphelie 
+      PRINT*,'obliquite en degres  :',obliquit
+      PRINT*,'Jours dans l annee : ',year_day
+      PRINT*,'Date perihelie : ',peri_day
+      unitastr=149.597870
+      e_elips=(aphelie-periheli)/(periheli+aphelie)
+      p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
+
+      print*,'e_elips',e_elips
+      print*,'p_elips',p_elips
+
+c-----------------------------------------------------------------------
+c calcul de l'angle polaire et de la distance au soleil :
+c -------------------------------------------------------
+
+c  calcul de l'zanomalie moyenne
+
+      zz=(year_day-pperi_day)/year_day
+      zanom=2.*pi*(zz-nint(zz))
+      zxref=abs(zanom)
+      PRINT*,'zanom  ',zanom
+
+c  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
+c  methode de Newton
+
+      zx0=zxref+R_ecc*sin(zxref)
+      DO 110 iter=1,100
+         zdx=-(zx0-R_ecc*sin(zx0)-zxref)/(1.-R_ecc*cos(zx0))
+         if(abs(zdx).le.(1.e-12)) goto 120
+         zx0=zx0+zdx
+110   continue
+120   continue
+      zx0=zx0+zdx
+      if(zanom.lt.0.) zx0=-zx0
+      PRINT*,'zx0   ',zx0
+
+c zteta est la longitude solaire
+
+      timeperi=2.*atan(sqrt((1.+R_ecc)/(1.-R_ecc))*tan(zx0/2.))
+      PRINT*,'longitude solaire du perihelie timeperi = ',timeperi
+
+      RETURN
+      END
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1198)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1201)
@@ -4,5 +4,5 @@
 
       SUBROUTINE physiq (nlon,nlev,
-     .            debut,lafin,rjourvrai,gmtime,pdtphys,
+     .            debut,lafin,jD_cur, jH_cur,pdtphys,
      .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
      .            u,v,t,qx,
@@ -60,6 +60,6 @@
 c debut---input-L-variable logique indiquant le premier passage
 c lafin---input-L-variable logique indiquant le dernier passage
-c rjour---input-R-numero du jour de l'experience
-c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
+c jD_cur--input-R-jour courant a l'appel de la physique (jour julien)
+c jH_cur--input-R-heure courante a l'appel de la physique (jour julien)
 c pdtphys-input-R-pas d'integration pour la physique (seconde)
 c paprs---input-R-pression pour chaque inter-couche (en Pa)
@@ -186,6 +186,6 @@
       INTEGER nlon
       INTEGER nlev
-      REAL rjourvrai
-      REAL gmtime
+      REAL :: jD_cur, jH_cur
+
       REAL pdtphys
       LOGICAL debut, lafin
@@ -703,5 +703,16 @@
 c Conditions aux limites
 c
-      INTEGER julien
+!
+! Gestion calendrier
+!
+      REAL :: jD_1jan, jH_1jan
+      INTEGER :: year_cur, mth_cur, day_cur, days_elapsed
+      REAL :: hour, day_since_equinox
+! Date de l'equinoxe de printemps
+      INTEGER, parameter :: mth_eq=3, day_eq=21
+      REAL :: jD_eq
+
+      LOGICAL, parameter :: new_orbit = .true.
+
 c
       INTEGER lmt_pas
@@ -1087,4 +1098,6 @@
 c$OMP THREADPRIVATE(first)
 
+      integer iunit
+
       logical, save::  read_climoz ! read ozone climatology
       integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
@@ -1114,7 +1127,7 @@
          write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
          write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
+     s 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
          write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys 
+     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys 
 
          write(lunout,*) 'papers, play, phi, u, v, t, omega'
@@ -1184,5 +1197,13 @@
 
 c======================================================================
-      xjour = rjourvrai
+! Gestion calendrier
+!
+      call ju2ymds(jD_cur+jH_cur, year_cur, mth_cur, day_cur, hour)
+      call ymds2ju(year_cur, 1, 1, 0., jD_1jan)
+      jH_1jan = jD_1jan - int (jD_1jan)
+      jD_1jan = int (jD_1jan) 
+      xjour = jD_cur - jD_1jan
+      days_elapsed = jD_cur - jD_1jan
+
 c
 c Si c'est le debut, il faut initialiser plusieurs choses
@@ -1472,6 +1493,5 @@
 
 cXXXPB Positionner date0 pour initialisation de ORCHIDEE
-      date0 = zjulian
-C      date0 = day_ini
+      date0 = jD_ref 
       WRITE(*,*) 'physiq date0 : ',date0
 c
@@ -1491,7 +1511,8 @@
          CALL VTe(VTphysiq)
          CALL VTb(VTinca)
-         iii = MOD(NINT(xjour),360)
-         calday = FLOAT(iii) + gmtime
-         WRITE(lunout,*) 'initial time ', xjour, calday
+!         iii = MOD(NINT(xjour),360)
+!         calday = FLOAT(iii) + jH_cur
+         calday = FLOAT(days_elapsed) + jH_cur
+         WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
 
          CALL chemini( 
@@ -1534,7 +1555,4 @@
 !
       itap   = itap + 1
-      julien = MOD(NINT(xjour),360)
-      if (julien .eq. 0) julien = 360
-
 !
 ! Update fraction of the sub-surfaces (pctsrf) and 
@@ -1542,7 +1560,6 @@
 ! on the surface fraction.
 !
-      CALL change_srf_frac(itap, dtime, julien, 
+      CALL change_srf_frac(itap, dtime, days_elapsed+1, 
      *     pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
-
 
 ! Tendances bidons pour les processus qui n'affectent pas certaines
@@ -1697,5 +1714,6 @@
          if (read_climoz) then
 C           Ozone climatology from a NetCDF file
-            call regr_pr_av(ncid_climoz, "tro3", julien, press_climoz,
+            call regr_pr_av(ncid_climoz, "tro3", days_elapsed+1, 
+     &           press_climoz,
      $           paprs, wo)
 !           Convert from mole fraction of ozone to column density of ozone in a
@@ -1707,5 +1725,6 @@
 C           "zmasse" changes a little.)
          else
-            CALL ozonecm(real(julien), rlat, paprs, wo)
+            CALL ozonecm(real(days_elapsed+1), rlat, paprs, wo)
+
          end if
       ENDIF
@@ -1749,20 +1768,28 @@
 ! doit donc etre placé avant radlwsw et pbl_surface
 
+! calcul selon la routine utilisee pour les planetes
+      if (new_orbit) then
+        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
+        day_since_equinox = (jD_cur + jH_cur) - jD_eq 
+!        day_since_equinox = (jD_cur) - jD_eq 
+        call solarlong(day_since_equinox, zlongi, dist)
+      else      
+! calcul selon la routine utilisee pour l'AR4
 !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
 !   solarlong0
-
-      if (solarlong0<-999.) then
-         CALL orbite(FLOAT(julien),zlongi,dist)
-      else
-         zlongi=solarlong0  ! longitude solaire vraie
-         dist=1.            ! distance au soleil / moyenne
+        if (solarlong0<-999.) then
+           CALL orbite(FLOAT(days_elapsed+1),zlongi,dist)
+        else
+           zlongi=solarlong0  ! longitude solaire vraie
+           dist=1.            ! distance au soleil / moyenne
+        endif
       endif
-
-      if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
+      if(prt_level.ge.1)                                                &
+     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
 
 !  Avec ou sans cycle diurne
       IF (cycle_diurne) THEN
         zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
-        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
+        CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
       ELSE
         CALL angle(zlongi, rlat, fract, rmu0)
@@ -1797,5 +1824,5 @@
 
       CALL pbl_surface( 
-     e     dtime,     date0,     itap,    julien,
+     e     dtime,     date0,     itap,    days_elapsed+1,
      e     debut,     lafin,
      e     rlon,      rlat,      rugoro,  rmu0,     
@@ -2094,6 +2121,8 @@
 
             if (itop_con(i).gt.klev-3) then
-               print*,'La convection monte trop haut '
-               print*,'itop_con(,',i,',)=',itop_con(i)
+              if(prt_level >= 9) then
+                write(lunout,*)'La convection monte trop haut '
+                write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
+              endif
             endif
           ENDDO
@@ -2635,5 +2664,5 @@
          IF (.NOT. aerosol_couple)
      &        CALL readaerosol_optic(
-     &        debut, new_aod, flag_aerosol, rjourvrai, pdtphys,
+     &        debut, new_aod, flag_aerosol, jD_cur-jD_ref, pdtphys,
      &        pplay, paprs, t_seri, rhcl, 
      &        mass_solu_aero, mass_solu_aero_pi, 
@@ -2764,5 +2793,5 @@
          CALL VTe(VTphysiq)
          CALL VTb(VTinca)
-         calday = FLOAT(julien) + gmtime
+         calday = FLOAT(days_elapsed + 1) + jH_cur
 
          IF (config_inca == 'aero') THEN
@@ -2775,6 +2804,6 @@
 
          CALL chemhook_begin (calday,
-     $                          julien,
-     $                          gmtime,
+     $                          days_elapsed,
+     $                          jH_cur,
      $                          pctsrf(1,1),
      $                          rlat,
@@ -3106,5 +3135,5 @@
       IF (is_sequential) THEN
       
-        CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
+        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,
      C                 ra,rg,romega,
      C                 rlat,rlon,pphis,
@@ -3133,5 +3162,5 @@
 
       call phytrac (
-     I     itap,     julien,    gmtime,   debut,
+     I     itap,     days_elapsed+1,    jH_cur,   debut,
      I     lafin,    dtime,     u, v,     t,
      I     paprs,    pplay,     pmfu,     pmfd, 
@@ -3351,7 +3380,7 @@
       write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
       write(lunout,*)
-     s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
+     s 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
       write(lunout,*)
-     s  nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
+     s  nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys,
      s  pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
      s  pctsrf(igout,is_sic)
@@ -3447,4 +3476,5 @@
       ENDIF
       
+!      first=.false.
 
       RETURN
Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/solarlong.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/solarlong.F	(revision 1201)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/solarlong.F	(revision 1201)
@@ -0,0 +1,131 @@
+      SUBROUTINE solarlong(pday,psollong,pdist_sol)
+
+      USE ioipsl
+
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Objet:
+c   ------
+c
+c      Calcul de la distance soleil-planete et de la declinaison
+c   en fonction du jour de l'annee.
+c
+c
+c   Methode:
+c   --------
+c
+c      Calcul complet de l'elipse
+c
+c   Interface:
+c   ----------
+c
+c      Uncommon comprenant les parametres orbitaux.
+c
+c   Arguments:
+c   ----------
+c
+c   Input:
+c   ------
+c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
+c   lwrite        clef logique pour sorties de controle
+c
+c   Output:
+c   -------
+c   pdist_sol     distance entre le soleil et la planete
+c                 ( en unite astronomique pour utiliser la constante 
+c                  solaire terrestre 1370 Wm-2 )
+c   pdecli        declinaison ( en radians )
+c
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+#include "planete.h"
+#include "YOMCST.h"
+      include 'iniprint.h'
+
+c arguments:
+c ----------
+
+      REAL pday,pdist_sol,pdecli,psollong
+      LOGICAL lwrite
+
+c Local:
+c ------
+
+      REAL zanom,xref,zx0,zdx,zteta,zz,pi
+      INTEGER iter
+      REAL :: pyear_day,pperi_day
+      REAL :: jD_eq, jD_peri
+
+c-----------------------------------------------------------------------
+c calcul de l'angle polaire et de la distance au soleil :
+c -------------------------------------------------------
+
+c   Initialisation eventuelle:
+      if(.not.unitastr.gt.1.e-4) then
+        call ioget_calendar(pyear_day)
+        call ymds2ju(2000, 3, 21, 0., jD_eq)
+        call ymds2ju(2001, 1, 4, 0., jD_peri)
+        pperi_day = jD_peri - jD_eq
+        pperi_day = R_peri + 180.
+        write(lunout,*)' Number of days in a year = ',pyear_day
+c         call iniorbit(249.22,206.66,669.,485.,25.2)
+         call iniorbit(152.59,146.61,pyear_day,pperi_day,R_incl)
+      endif
+
+c  calcul de l'zanomalie moyenne
+
+      zz=(pday-peri_day)/year_day
+      pi=2.*asin(1.)
+      zanom=2.*pi*(zz-nint(zz))
+      xref=abs(zanom)
+
+c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
+c  methode de Newton
+
+!      zx0=xref+e_elips*sin(xref)
+      zx0=xref+R_ecc*sin(xref)
+      DO 110 iter=1,10
+!         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
+         zdx=-(zx0-R_ecc*sin(zx0)-xref)/(1.-R_ecc*cos(zx0))
+         if(abs(zdx).le.(1.e-7)) goto 120
+         zx0=zx0+zdx
+110   continue
+120   continue
+      zx0=zx0+zdx
+      if(zanom.lt.0.) zx0=-zx0
+
+c zteta est la longitude solaire
+
+!      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+      zteta=2.*atan(sqrt((1.+R_ecc)/(1.-R_ecc))*tan(zx0/2.))
+
+      psollong=zteta-timeperi
+
+      IF(psollong.LT.0.) psollong=psollong+2.*pi
+      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
+
+      psollong = psollong * 180. / pi
+
+c distance soleil
+
+      pdist_sol = (1-R_ecc*R_ecc)
+     &      /(1+R_ecc*COS(pi/180.*(psollong-(R_peri+180.0))))
+!      pdist_sol = (1-e_elips*e_elips)
+!     &      /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0))))
+c-----------------------------------------------------------------------
+c   sorties eventuelles:
+c   ---------------------
+
+c     IF (lwrite) THEN
+c        PRINT*,'jour de l"annee   :',pday
+c        PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
+c        PRINT*,'declinaison (en degres) :',pdecli*180./pi
+c     ENDIF
+
+      RETURN
+      END
