Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iostart.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iostart.F90	(revision 269)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iostart.F90	(revision 270)
@@ -44,9 +44,10 @@
 CONTAINS
 
-  SUBROUTINE open_startphy(filename)
+  SUBROUTINE open_startphy(filename,found)
   USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
   USE mod_phys_lmdz_para, only: is_master, bcast
   IMPLICIT NONE
-    CHARACTER(LEN=*) :: filename
+    CHARACTER(LEN=*),INTENT(IN) :: filename
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
     INTEGER          :: ierr
 
@@ -56,9 +57,13 @@
         write(*,*)'open_startphy: problem opening file '//trim(filename)
         write(*,*)trim(nf90_strerror(ierr))
-        CALL ABORT
+        !CALL ABORT
+        found=.false.
+      ELSE
+        found=.true.
       ENDIF
     ENDIF
     
     CALL bcast(nid_start) ! tell all procs about nid_start
+    if(present(found)) CALL bcast(found)
   
   END SUBROUTINE open_startphy
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90	(revision 269)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90	(revision 270)
@@ -83,5 +83,5 @@
       
       INTEGER :: indextime=1 ! index of selected time, default value=1
-      logical :: found
+      logical :: found,found_file
 
 !
@@ -98,13 +98,18 @@
 
 ! open physics initial state file:
-call open_startphy(fichnom)
-
-
-! possibility to modify tab_cntrl in tabfi
-write(*,*)
-write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
-call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
+call open_startphy(fichnom,found_file)
+
+! Ehouarn, if file not found, then call tabfi with nid_start==0
+if (.not.found_file) then
+  write(*,*) 'phyetat0_academic: call tabfi with nid_start=0'
+  call tabfi (ngrid,0,Lmodif,tab0,day_ini,lmax,p_rad, &
                    p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
-
+else
+  ! possibility to modify tab_cntrl in tabfi
+  write(*,*)
+  write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
+  call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
+                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+endif
 !c
 !c Lecture des latitudes (coordonnees):
@@ -165,5 +170,9 @@
 
 ! Load surface geopotential:
-call get_field("phisfi",phisfi,found)
+if (found_file) then
+  call get_field("phisfi",phisfi,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <phisfi>"
@@ -175,5 +184,9 @@
 
 ! Load bare ground albedo:
-call get_field("albedodat",albedodat,found)
+if (found_file) then
+  call get_field("albedodat",albedodat,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <albedodat>"
@@ -187,5 +200,9 @@
 
 ! ZMEA
-call get_field("ZMEA",zmea,found)
+if (found_file) then
+  call get_field("ZMEA",zmea,found)
+else
+  found=.false.
+endif
 if (.not.found) then
     zmea(:)=0.
@@ -196,5 +213,9 @@
 
 ! ZSTD
-call get_field("ZSTD",zstd,found)
+if (found_file) then
+  call get_field("ZSTD",zstd,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <ZSTD>"
@@ -206,5 +227,9 @@
 
 ! ZSIG
-call get_field("ZSIG",zsig,found)
+if (found_file) then
+  call get_field("ZSIG",zsig,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <ZSIG>"
@@ -216,5 +241,9 @@
 
 ! ZGAM
-call get_field("ZGAM",zgam,found)
+if (found_file) then
+  call get_field("ZGAM",zgam,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <ZGAM>"
@@ -226,5 +255,9 @@
 
 ! ZTHE
-call get_field("ZTHE",zthe,found)
+if (found_file) then
+  call get_field("ZTHE",zthe,found)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <ZTHE>"
@@ -236,9 +269,12 @@
 
 ! Surface temperature :
-call get_field("tsurf",tsurf,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tsurf>"
-!mi initialising tsurf with pt(:,1) 
-!tsurf(:)=175.0
+if (found_file) then
+  call get_field("tsurf",tsurf,found,indextime)
+else
+  found=.false.
+endif
+if (.not.found) then
+  !mi initialising tsurf with pt(:,1)
+  !tsurf(:)=175.0
 else
   write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
@@ -247,5 +283,9 @@
 
 ! Surface emissivity
-call get_field("emis",emis,found,indextime)
+if (found_file) then
+  call get_field("emis",emis,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <emis>"
@@ -257,5 +297,9 @@
 
 ! Cloud fraction (added by BC 2010)
-call get_field("cloudfrac",cloudfrac,found,indextime)
+if (found_file) then
+  call get_field("cloudfrac",cloudfrac,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <cloudfrac>"
@@ -267,5 +311,9 @@
 
 ! Total cloud fraction (added by BC 2010)
-call get_field("totcloudfrac",totcloudfrac,found,indextime)
+if (found_file) then
+  call get_field("totcloudfrac",totcloudfrac,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <totcloudfrac>"
@@ -277,5 +325,9 @@
 
 ! Height of oceanic ice (added by BC 2010)
-call get_field("hice",hice,found,indextime)
+if (found_file) then
+  call get_field("hice",hice,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <hice>"
@@ -291,5 +343,9 @@
 ! SLAB OCEAN (added by BC 2014)
 ! nature of the surface
-call get_field("rnat",rnat,found,indextime)
+if (found_file) then
+  call get_field("rnat",rnat,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <rnat>"
@@ -310,5 +366,9 @@
 endif
 ! Pourcentage of sea ice cover
-call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
+if (found_file) then
+  call get_field("pctsrf_sic",pctsrf_sic,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
@@ -321,5 +381,9 @@
 endif
 ! Slab ocean temperature (2 layers)
-call get_field("tslab",tslab,found,indextime)
+if (found_file) then
+  call get_field("tslab",tslab,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <tslab>"
@@ -334,5 +398,9 @@
 endif
 ! Oceanic ice temperature
-call get_field("tsea_ice",tsea_ice,found,indextime)
+if (found_file) then
+  call get_field("tsea_ice",tsea_ice,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <tsea_ice>"
@@ -345,5 +413,9 @@
 endif
 !  Oceanic ice quantity (kg/m^2)
-call get_field("sea_ice",sea_ice,found,indextime)
+if (found_file) then
+  call get_field("sea_ice",sea_ice,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <sea_ice>"
@@ -360,5 +432,9 @@
 
 ! pbl wind variance
-call get_field("q2",q2,found,indextime)
+if (found_file) then
+  call get_field("q2",q2,found,indextime)
+else
+  found=.false.
+endif
 if (.not.found) then
   write(*,*) "phyetat0: Failed loading <q2>"
@@ -380,5 +456,9 @@
                            ' h2o_ice instead of h2o_vap'
     endif
-    call get_field(txt,qsurf(:,iq),found,indextime)
+    if (found_file) then
+      call get_field(txt,qsurf(:,iq),found,indextime)
+    else
+      found=.false.
+    endif
     if (.not.found) then
       write(*,*) "phyetat0: Failed loading <",trim(txt),">"
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 269)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 270)
@@ -569,4 +569,5 @@
 !ym               rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
+         write(*,*) "physiq: firstcall, call phyetat0_academic"
          call phyetat0_academic(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,   &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
@@ -575,4 +576,5 @@
 
 !mi initialising tsurf with pt
+         write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
          tsurf(:)=pt(:,1)
 
@@ -744,7 +746,8 @@
 
          if (ngrid.ne.1) then ! no need to create a restart file in 1d
-           call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
-                         ptimestep,pday+nday,time_phys,area, &
-                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+! EM: No restart file (for now).
+!           call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
+!                         ptimestep,pday+nday,time_phys,area, &
+!                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
          endif
          
@@ -1972,9 +1975,11 @@
 !                    cloudfrac,totcloudfrac,hice,noms)
 !#endif
-              call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
-                      ptimestep,ztime_fin, &
-                      tsurf,tsoil,emis,q2,qsurf_hist, &
-                      cloudfrac,totcloudfrac,hice, &
-                      rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
+
+! EM: do not write a restart file (for now).
+!              call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
+!                      ptimestep,ztime_fin, &
+!                      tsurf,tsoil,emis,q2,qsurf_hist, &
+!                      cloudfrac,totcloudfrac,hice, &
+!                      rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
             endif
 
@@ -2239,4 +2244,14 @@
       CALL write_xios_field("ASR",fluxabs_sw)
       CALL write_xios_field("OLR",fluxtop_lw)
+      call write_xios_field("input_temp",pt)
+      call write_xios_field("input_u",pu)
+      call write_xios_field("input_v",pv)
+      call write_xios_field("dtrad",dtrad)
+      call write_xios_field("zdtlw",zdtlw)
+      call write_xios_field("zdtsw",zdtsw)
+      call write_xios_field("zdtdyn",zdtdyn/ptimestep)
+      call write_xios_field("zdtdif",zdtdif)
+      call write_xios_field("zdtadj",zdtadj)
+      call write_xios_field("pdt",pdt)
       IF (lastcall) CALL finalize_xios_output
       
@@ -2344,4 +2359,5 @@
       endif
 
+      write(*,*) "physiq: done, zday=",zday
       return
     end subroutine physiq
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F	(revision 269)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F	(revision 270)
@@ -44,5 +44,5 @@
 c=======================================================================
 ! to use  'getin'
-      use ioipsl_getincom , only: getin
+      use ioipsl_getincom_p , only: getin_p
 
       use surfdat_h, only: albedice, emisice, iceradius, dtemisice,
@@ -88,4 +88,61 @@
 c  Initialization of various physical constants to defaut values (nid = 0 case)
 c-----------------------------------------------------------------------
+        ! Ehouarn: Default Saturn values:
+        tab_cntrl(:)=0
+        write(*,*) "Using default Saturn values..."
+        ! these should be read in a def file I guess...
+        lmax=0 ! not used anyways
+        !day_ini=0
+        time=0
+        ! radius of the planet
+        rad=60268000
+        call getin_p("radius",rad)
+        ! Planetary rotation rate
+        omeg=0.00016512100410182
+        call getin_p("omega",omeg)
+        ! Gravity
+        g=10.44
+        call getin_p("g",g)
+        !mugaz=2.34 !EM: does not give cpp=11500
+        mugaz=2.53 ! with this value of mugaz, cpp=11500
+        call getin_p("mugaz",mugaz)
+        ! kappa
+        rcp=0.2857143
+        call getin_p("kappa",rcp)
+        cpp=(8.314511/(mugaz/1000.0))/rcp
+        call getin_p("cpp",cpp)
+!        write(*,*) "tabfi: cpp=",cpp
+        ! length (s) of a "standard" day
+        daysec=38052
+        call getin_p("day_length",daysec)
+        ! physics time step (s) ! not sure we need this here
+        dtphys=19026
+        ! length of year, in standard days
+        year_day=24430
+        ! Orbital parameters
+        periastr=9.02151966094971
+        apoastr=10.054479598999
+        peri_day=19280
+        obliquit=26.7299995422363
+        ! Other parameters some physical parÃ©ametrizations need
+        z0=1e-2
+        lmixmin=30
+        emin_turb=1.e-6
+        albedice(:)=0
+        emisice(:)=0
+        emissiv=0
+        iceradius(:)=1.e-6
+        dtemisice(:)=0
+        volcapa=1000000
+c-----------------------------------------------------------------------
+c	Save some constants for later use (as routine arguments)
+c-----------------------------------------------------------------------
+        p_omeg = omeg
+        p_g = g
+        p_cpp = cpp
+        p_mugaz = mugaz
+        p_daysec = daysec
+        p_rad=rad
+
       ELSE
 c-----------------------------------------------------------------------
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/xios_output.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/xios_output.F90	(revision 269)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/xios_output.F90	(revision 270)
@@ -20,4 +20,5 @@
   USE mod_phys_lmdz_para
   USE mod_grid_phy_lmdz
+  USE planete_mod, ONLY: ap,bp,preff
   
   IMPLICIT NONE
@@ -33,4 +34,5 @@
     REAL :: bounds_lat_mpi(klon_mpi,nvertex)
     INTEGER :: l
+    REAL :: presnivs(klev)
 
     CALL gather_omp(rlond,lon_mpi)
@@ -45,8 +47,12 @@
     CALL xios_set_current_context(ctx_hdl)
    
-    lev_value(:) = (/ (l,l=1,klev) /)     
+!    lev_value(:) = (/ (l,l=1,klev) /)     
+!    CALL xios_set_axis_attr("altitude",size=klev ,value=lev_value)
     
-    CALL xios_set_axis_attr("altitude",size=klev ,value=lev_value)
-    
+    ! build presnivs "pseudo-pressure" at mid-layers
+    DO l=1,klev
+      presnivs(l)=0.5*(ap(l)+bp(l)*preff+ap(l+1)+bp(l+1)*preff)
+    ENDDO
+    CALL xios_set_axis_attr("presnivs",size=klev,value=presnivs,unit="Pa",positive="down")
   
     CALL xios_set_domain_attr("horiz_domain",ni_glo=klon_glo, ibegin=klon_mpi_begin, ni=klon_mpi)
