Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1668)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1670)
@@ -35,5 +35,7 @@
       USE temps_mod, ONLY: day_ini
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      use tabfi_mod, only: tabfi
       use iniphysiq_mod, only: iniphysiq
+      use phyetat0_mod, only: phyetat0
       implicit none
 
@@ -337,5 +339,5 @@
         write(*,*) 'Reading file STARTFI'
         fichnom = 'startfi.nc'
-        CALL phyetat0 (ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
+        CALL phyetat0(.true.,ngridmx,llm,fichnom,tab0,Lmodif,nsoilmx,
      .        nqtot,day_ini,time,
      .        tsurf,tsoil,emis,q2,qsurf)
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 1668)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F	(revision 1670)
@@ -35,4 +35,5 @@
       USE temps_mod, ONLY: day_ini
       USE iniphysiq_mod, ONLY: iniphysiq
+      use phyetat0_mod, only: phyetat0
       implicit none
 
@@ -214,5 +215,5 @@
 
 
-      CALL phyetat0 (ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
+      CALL phyetat0(.true.,ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
      .      day_ini_fi,timefi,
      .      tsurf,tsoil,emis,q2,qsurf)
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1668)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1670)
@@ -74,3 +74,7 @@
 !$OMP THREADPRIVATE(iscallphys)
 
+      ! do we read a startphy.nc file (default=.true.)
+      logical,save :: startphy_file=.true. 
+!$OMP THREADPRIVATE(startphy_file)
+
 END MODULE callkeys_mod
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1668)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1670)
@@ -95,4 +95,7 @@
   call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
 
+  ! do we read a startphy.nc file? (default: .true.)
+  call getin_p("startphy_file",startphy_file)
+  
 ! --------------------------------------------------------------
 !  Reading the "callphys.def" file controlling some key options
Index: trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90	(revision 1668)
+++ 	(revision )
@@ -1,289 +1,0 @@
-subroutine phyetat0 (ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
-                     day_ini,time,tsurf,tsoil, &
-                     emis,q2,qsurf)
-
-
-  USE tracer_h, ONLY: noms
-  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
-  use iostart, only: nid_start, open_startphy, close_startphy, &
-                     get_field, get_var, inquire_field, &
-                     inquire_dimension, inquire_dimension_length
-
-  implicit none
-
-!======================================================================
-! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
-!  Adaptation à Mars : Yann Wanherdrick 
-! Objet: Lecture de l etat initial pour la physique
-!======================================================================
-#include "netcdf.inc"
-
-!======================================================================
-!  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
-!  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
-!======================================================================
-!  Arguments:
-!  ---------
-!  inputs:
-  integer,intent(in) :: ngrid
-  integer,intent(in) :: nlayer
-  character*(*),intent(in) :: fichnom ! "startfi.nc" file
-  integer,intent(in) :: tab0
-  integer,intent(in) :: Lmodif
-  integer,intent(in) :: nsoil ! # of soil layers
-  integer,intent(in) :: nq
-  integer,intent(in) :: day_ini
-  real,intent(in) :: time
-
-!  outputs:
-  real,intent(out) :: tsurf(ngrid) ! surface temperature
-  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
-  real,intent(out) :: emis(ngrid) ! surface emissivity
-  real,intent(out) :: q2(ngrid,nlayer+1) ! 
-  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
-
-!======================================================================
-!  Local variables:
-
-!      INTEGER radpas
-!      REAL solaire
-
-      real xmin,xmax ! to display min and max of a field
-!
-      INTEGER ig,iq,lmax
-      INTEGER nid, nvarid
-      INTEGER ierr, i, nsrf
-!      integer isoil 
-!      INTEGER length
-!      PARAMETER (length=100)
-      CHARACTER*7 str7
-      CHARACTER*2 str2
-      CHARACTER*1 yes
-!
-      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
-      INTEGER nqold
-
-! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
-!      logical :: oldtracernames=.false.
-      integer :: count
-      character(len=30) :: txt ! to store some text
-      
-      INTEGER :: indextime=1 ! index of selected time, default value=1
-      logical :: found
-
-!
-! ALLOCATE ARRAYS IN surfdat_h
-!
-IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
-IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
-IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
-IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
-IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
-IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
-IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
-
-
-! 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, &
-                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
-
-!c
-!c Lecture des latitudes (coordonnees):
-!c
-!      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
-!         CALL abort
-!      ENDIF
-!#ifdef NC_DOUBLE
-!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
-!#else
-!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
-!#endif
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
-!         CALL abort
-!      ENDIF
-!c
-!c Lecture des longitudes (coordonnees):
-!c
-!      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
-!         CALL abort
-!      ENDIF
-!#ifdef NC_DOUBLE
-!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
-!#else
-!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
-!#endif
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
-!         CALL abort
-!      ENDIF
-!c
-!c Lecture des aires des mailles:
-!c
-!      ierr = NF_INQ_VARID (nid, "area", nvarid)
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Le champ <area> est absent'
-!         CALL abort
-!      ENDIF
-!#ifdef NC_DOUBLE
-!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
-!#else
-!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
-!#endif
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
-!         CALL abort
-!      ENDIF
-!      xmin = 1.0E+20
-!      xmax = -1.0E+20
-!      xmin = MINVAL(area)
-!      xmax = MAXVAL(area)
-!      PRINT*,'Aires des mailles <area>:', xmin, xmax
-
-! Load surface geopotential:
-call get_field("phisfi",phisfi,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <phisfi>"
-  call abort
-else
-  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
-             minval(phisfi), maxval(phisfi)
-endif
-
-! Load bare ground albedo:
-call get_field("albedodat",albedodat,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <albedodat>"
-  call abort
-else
-  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
-             minval(albedodat), maxval(albedodat)
-endif
-
-! ZMEA
-call get_field("ZMEA",zmea,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZMEA>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZMEA> range:", &
-             minval(zmea), maxval(zmea)
-endif
-
-! ZSTD
-call get_field("ZSTD",zstd,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSTD>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZSTD> range:", &
-             minval(zstd), maxval(zstd)
-endif
-
-! ZSIG
-call get_field("ZSIG",zsig,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSIG>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZSIG> range:", &
-             minval(zsig), maxval(zsig)
-endif
-
-! ZGAM
-call get_field("ZGAM",zgam,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZGAM>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZGAM> range:", &
-             minval(zgam), maxval(zgam)
-endif
-
-! ZTHE
-call get_field("ZTHE",zthe,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZTHE>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZTHE> range:", &
-             minval(zthe), maxval(zthe)
-endif
-
-! Surface temperature :
-call get_field("tsurf",tsurf,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tsurf>"
-  call abort
-else
-  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
-             minval(tsurf), maxval(tsurf)
-endif
-
-! Surface emissivity
-call get_field("emis",emis,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <emis>"
-  call abort
-else
-  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
-             minval(emis), maxval(emis)
-endif
-
-! pbl wind variance
-call get_field("q2",q2,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <q2>"
-  call abort
-else
-  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
-             minval(q2), maxval(q2)
-endif
-
-! tracer on surface
-if (nq.ge.1) then
-  do iq=1,nq
-    txt=noms(iq)
-    
-    !! There was a bug here. MT2015.
-    
-    !if (txt.eq."h2o_vap") then
-      ! There is no surface tracer for h2o_vap;
-      ! "h2o_ice" should be loaded instead
-     ! txt="h2o_ice"
-     ! write(*,*) 'phyetat0: loading surface tracer', &
-     !                      ' h2o_ice instead of h2o_vap'
-    !endif
-    
-    call get_field(txt,qsurf(:,iq),found,indextime)
-    if (.not.found) then
-      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
-      write(*,*) "         ",trim(txt)," is set to zero"
-      qsurf(:,iq) = 0.
-    else
-      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
-                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
-    endif
-  enddo
-endif ! of if (nq.ge.1)
-
-
-! Call to soil_settings, in order to read soil temperatures,
-! as well as thermal inertia and volumetric heat capacity
-call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
-!
-! close file:
-!
-call close_startphy
-
-END SUBROUTINE phyetat0
Index: trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 1670)
+++ trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 1670)
@@ -0,0 +1,253 @@
+module phyetat0_mod
+
+implicit none
+
+contains
+
+subroutine phyetat0 (startphy_file, &
+                     ngrid,nlayer,fichnom,tab0,Lmodif,nsoil,nq, &
+                     day_ini,time,tsurf,tsoil, &
+                     emis,q2,qsurf)
+
+
+  use tabfi_mod, only: tabfi
+  USE tracer_h, ONLY: noms
+  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
+  use iostart, only: nid_start, open_startphy, close_startphy, &
+                     get_field, get_var, inquire_field, &
+                     inquire_dimension, inquire_dimension_length
+
+  implicit none
+
+!======================================================================
+!  Arguments:
+!  ---------
+!  inputs:
+  logical,intent(in) :: startphy_file ! .true. if reading start file
+  integer,intent(in) :: ngrid
+  integer,intent(in) :: nlayer
+  character*(*),intent(in) :: fichnom ! "startfi.nc" file
+  integer,intent(in) :: tab0
+  integer,intent(in) :: Lmodif
+  integer,intent(in) :: nsoil ! # of soil layers
+  integer,intent(in) :: nq
+  integer,intent(out) :: day_ini
+  real,intent(out) :: time
+
+!  outputs:
+  real,intent(out) :: tsurf(ngrid) ! surface temperature
+  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
+  real,intent(out) :: emis(ngrid) ! surface emissivity
+  real,intent(out) :: q2(ngrid,nlayer+1) ! 
+  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+
+!======================================================================
+!  Local variables:
+
+!      INTEGER radpas
+!      REAL solaire
+
+      real xmin,xmax ! to display min and max of a field
+!
+      INTEGER ig,iq,lmax
+      INTEGER nid, nvarid
+      INTEGER ierr, i, nsrf
+!      integer isoil 
+!      INTEGER length
+!      PARAMETER (length=100)
+      CHARACTER*7 str7
+      CHARACTER*2 str2
+      CHARACTER*1 yes
+!
+      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
+      INTEGER nqold
+
+! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
+!      logical :: oldtracernames=.false.
+      integer :: count
+      character(len=30) :: txt ! to store some text
+      
+      INTEGER :: indextime=1 ! index of selected time, default value=1
+      logical :: found
+      
+      character(len=8) :: modname="phyetat0"
+
+!
+! ALLOCATE ARRAYS IN surfdat_h
+!
+IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
+IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
+IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
+IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
+IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
+IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
+IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
+
+if (startphy_file) then
+  ! 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, &
+                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+
+else ! "academic" initialization of planetary parameters via tabfi
+  call tabfi (ngrid,0,0,0,day_ini,lmax,p_rad, &
+                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+endif ! of if (startphy_file)
+
+if (startphy_file) then
+  ! Load surface geopotential:
+  call get_field("phisfi",phisfi,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <phisfi>",1)
+  endif
+else
+  phisfi(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
+               minval(phisfi), maxval(phisfi)
+
+if (startphy_file) then
+  ! Load bare ground albedo:
+  call get_field("albedodat",albedodat,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <albedodat>",1)
+  endif
+else
+  albedodat(:)=0.5 ! would be better to read value from def file...
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
+             minval(albedodat), maxval(albedodat)
+
+! ZMEA
+if (startphy_file) then
+  call get_field("ZMEA",zmea,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <ZMEA>",1)
+  endif
+else
+  zmea(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: <ZMEA> range:", &
+             minval(zmea), maxval(zmea)
+
+! ZSTD
+if (startphy_file) then
+  call get_field("ZSTD",zstd,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <ZSTD>",1)
+  endif
+else
+  zstd(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: <ZSTD> range:", &
+             minval(zstd), maxval(zstd)
+
+! ZSIG
+if (startphy_file) then
+  call get_field("ZSIG",zsig,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <ZSIG>",1)
+  endif
+else
+  zsig(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: <ZSIG> range:", &
+             minval(zsig), maxval(zsig)
+
+! ZGAM
+if (startphy_file) then
+  call get_field("ZGAM",zgam,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <ZGAM>",1)
+  endif
+else
+  zgam(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: <ZGAM> range:", &
+             minval(zgam), maxval(zgam)
+
+! ZTHE
+if (startphy_file) then
+  call get_field("ZTHE",zthe,found)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <ZTHE>",1)
+  endif
+else
+  zthe(:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: <ZTHE> range:", &
+             minval(zthe), maxval(zthe)
+
+! Surface temperature :
+if (startphy_file) then
+  call get_field("tsurf",tsurf,found,indextime)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <tsurf>",1)
+  endif
+else
+  tsurf(:)=0 ! will be updated afterwards in physiq !
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
+             minval(tsurf), maxval(tsurf)
+
+! Surface emissivity
+if (startphy_file) then
+  call get_field("emis",emis,found,indextime)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <emis>",1)
+  endif
+else
+  emis(:)=1 ! would be better to read value from def file...
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: Surface emissivity <emis> range:", &
+             minval(emis), maxval(emis)
+
+! pbl wind variance
+if (startphy_file) then
+  call get_field("q2",q2,found,indextime)
+  if (.not.found) then
+    call abort_physic(modname,"Failed loading <q2>",1)
+  endif
+else
+  q2(:,:)=0
+endif ! of if (startphy_file)
+write(*,*) "phyetat0: PBL wind variance <q2> range:", &
+             minval(q2), maxval(q2)
+
+! tracer on surface
+if (nq.ge.1) then
+  do iq=1,nq
+    txt=noms(iq)
+    if (startphy_file) then
+      call get_field(txt,qsurf(:,iq),found,indextime)
+      if (.not.found) then
+        write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+        write(*,*) "         ",trim(txt)," is set to zero"
+        qsurf(:,iq) = 0.
+      endif
+    else
+      qsurf(:,iq)=0
+    endif ! of if (startphy_file)
+    write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
+                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+  enddo! of do iq=1,nq
+endif ! of if (nq.ge.1)
+
+
+if (startphy_file) then
+  ! Call to soil_settings, in order to read soil temperatures,
+  ! as well as thermal inertia and volumetric heat capacity
+  call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+endif ! of if (startphy_file)
+!
+! close file:
+!
+if (startphy_file) call close_startphy
+
+end subroutine phyetat0
+
+end module phyetat0_mod
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1668)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1670)
@@ -25,4 +25,5 @@
                           alpha_lift, alpha_devil, qextrhor
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
+      use phyetat0_mod, only: phyetat0
       use phyredem, only: physdem0, physdem1
       use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
@@ -146,5 +147,5 @@
 !    ------------------
 
-#include "netcdf.inc"
+include "netcdf.inc"
 
 ! Arguments :
@@ -225,5 +226,5 @@
       real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
 
-      integer l,ig,ierr,iq,nw
+      integer l,ig,ierr,iq,nw,isoil
       
       ! FOR DIAGNOSTIC :
@@ -432,6 +433,18 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
+         call phyetat0(startphy_file,                                 &
+                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
                        day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
+         if (.not.startphy_file) then
+           ! additionnal "academic" initialization of physics
+           if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
+           tsurf(:)=pt(:,1)
+           if (is_master) write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
+           do isoil=1,nsoilmx
+             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
+           enddo
+           if (is_master) write(*,*) "Physiq: initializing day_ini to pdat !"
+           day_ini=pday
+         endif
 
          if (pday.ne.day_ini) then
Index: trunk/LMDZ.TITAN/libf/phytitan/tabfi.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tabfi.F	(revision 1668)
+++ 	(revision )
@@ -1,522 +1,0 @@
-c=======================================================================
-      SUBROUTINE tabfi(ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
-     .                 p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
-c=======================================================================
-c
-c   C. Hourdin 15/11/96
-c
-c   Object:        Lecture du tab_cntrl physique dans un fichier 
-c   ------            et initialisation des constantes physiques
-c
-c   Arguments:
-c   ----------
-c
-c     Inputs:
-c     ------
-c
-c      - nid:    unitne logique du fichier ou on va lire le tab_cntrl    
-c                      (ouvert dans le programme appellant) 
-c
-c                 si nid=0:
-c                       pas de lecture du tab_cntrl mais
-c                       Valeurs par default des constantes physiques
-c        
-c      - tab0:    Offset de tab_cntrl a partir duquel sont ranges 
-c                  les parametres physiques (50 pour start_archive)
-c
-c      - Lmodif:  si on souhaite modifier les constantes  Lmodif = 1 = TRUE
-c
-c
-c     Outputs:
-c     --------
-c
-c      - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite
-c                              comparer avec le day_ini dynamique)
-c
-c      - lmax:    tab_cntrl(tab0+2) (pour test avec nlayer)
-c
-c      - p_rad
-c      - p_omeg   !
-c      - p_g      ! Constantes physiques ayant des 
-c      - p_mugaz  ! homonymes dynamiques
-c      - p_daysec !
-c
-c=======================================================================
-! to use  'getin'
-      use ioipsl_getincom , only: getin
-
-      use surfdat_h, only: emisice, iceradius, dtemisice,
-     &                     emissiv
-      use comsoil_h, only: volcapa
-      use iostart, only: get_var
-      use mod_phys_lmdz_para, only: is_parallel
-      use planete_mod, only: year_day, periastr, apoastr, peri_day,
-     &                       obliquit, z0, lmixmin, emin_turb
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
-      use time_phylmdz_mod, only: dtphys, daysec
-      use callkeys_mod, only: check_cpp_match,force_cpp
-      implicit none
- 
-#include "netcdf.inc"
-
-c-----------------------------------------------------------------------
-c   Declarations
-c-----------------------------------------------------------------------
-
-c Arguments
-c ---------
-      INTEGER,INTENT(IN) :: ngrid,nid,tab0
-      INTEGER*4,INTENT(OUT) :: day_ini
-      INTEGER,INTENT(IN) :: Lmodif
-      INTEGER,INTENT(OUT) :: lmax
-      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
-
-c Variables
-c ---------
-      INTEGER,PARAMETER :: length=100
-      REAL tab_cntrl(length) ! array in which are stored the run's parameters
-      INTEGER  ierr,nvarid
-      INTEGER size
-      CHARACTER modif*20
-      LOGICAL :: found
-      
-      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
-
-      IF (nid.eq.0) then
-c-----------------------------------------------------------------------
-c  Initialization of various physical constants to defaut values (nid = 0 case)
-c-----------------------------------------------------------------------
-      ELSE
-c-----------------------------------------------------------------------
-c  Initialization of physical constants by reading array tab_cntrl(:)
-c		which contains these parameters	(nid != 0 case)
-c-----------------------------------------------------------------------
-c Read 'controle' array
-c
-
-       call get_var("controle",tab_cntrl,found)
-       if (.not.found) then
-         write(*,*)"tabfi: Failed reading <controle> array"
-         call abort
-       else
-         write(*,*)'tabfi: tab_cntrl',tab_cntrl
-       endif
-c
-c  Initialization of some physical constants
-c informations on physics grid
-!      if(ngrid.ne.tab_cntrl(tab0+1)) then
-!         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
-!         print*,tab_cntrl(tab0+1),ngrid
-!      endif
-      lmax = nint(tab_cntrl(tab0+2))
-      day_ini = tab_cntrl(tab0+3)
-      time = tab_cntrl(tab0+4)
-      write (*,*) 'IN tabfi day_ini=',day_ini
-c Informations about planet for dynamics and physics
-      rad = tab_cntrl(tab0+5)
-      omeg = tab_cntrl(tab0+6)
-      g = tab_cntrl(tab0+7)
-      mugaz = tab_cntrl(tab0+8)
-      rcp = tab_cntrl(tab0+9)
-      cpp=(8.314511/(mugaz/1000.0))/rcp
-      daysec = tab_cntrl(tab0+10)
-      dtphys = tab_cntrl(tab0+11)
-c Informations about planet for the physics only
-      year_day = tab_cntrl(tab0+14)
-      periastr = tab_cntrl(tab0+15)
-      apoastr = tab_cntrl(tab0+16)
-      peri_day = tab_cntrl(tab0+17)
-      obliquit = tab_cntrl(tab0+18)
-c boundary layer and turbeulence
-      z0 = tab_cntrl(tab0+19)
-      lmixmin = tab_cntrl(tab0+20)
-      emin_turb = tab_cntrl(tab0+21)
-c optical properties of polar caps and ground emissivity
-      emisice(1) = tab_cntrl(tab0+24)
-      emisice(2) = tab_cntrl(tab0+25)
-      emissiv    = tab_cntrl(tab0+26)
-      iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
-      iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
-      dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north)
-      dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south)
-c soil properties
-      volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
-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
-
-      ENDIF    ! end of (nid = 0) 
-
-c-----------------------------------------------------------------------
-c	Write physical constants to output before modifying them
-c-----------------------------------------------------------------------
- 
-   6  FORMAT(a20,e15.6,e15.6)
-   5  FORMAT(a20,f12.2,f12.2)
- 
-      write(*,*) '*****************************************************'
-      write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
-      write(*,*) '*****************************************************'
-      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
-      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
-      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
-      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
-      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
-      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
-      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
-      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
-      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
-      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
-
-      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
-      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
-      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
-      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
-      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
-
-      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
-      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
-      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
-
-      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
-      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
-      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
-      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
-      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
-      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
-      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
-
-      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
-
-      write(*,*)
-      write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif
-
-c-----------------------------------------------------------------------
-c	 Modifications...
-! NB: Modifying controls should only be done by newstart, and in seq mode
-      if ((Lmodif.eq.1).and.is_parallel) then
-        write(*,*) "tabfi: Error modifying tab_control should",
-     &             " only happen in serial mode (eg: by newstart)"
-        stop
-      endif
-c-----------------------------------------------------------------------
-
-      IF(Lmodif.eq.1) then
-
-      write(*,*)
-      write(*,*) 'Change values in tab_cntrl ? :'
-      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-      write(*,*) '(Current values given above)'
-      write(*,*)
-      write(*,*) '(3)          day_ini : Initial day (=0 at Ls=0)'
-      write(*,*) '(19)              z0 :  surface roughness (m)'
-      write(*,*) '(21)       emin_turb :  minimal energy (PBL)'
-      write(*,*) '(20)         lmixmin : mixing length (PBL)'
-      write(*,*) '(26)         emissiv : ground emissivity'
-      write(*,*) '(24 et 25)   emisice : CO2 ice max emissivity '
-      write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
-      write(*,*) '(33 et 34) dtemisice : time scale for snow',
-     &           'metamorphism'
-      write(*,*) '(35)      volcapa : soil volumetric heat capacity'
-      write(*,*) '(18)     obliquit : planet obliquity (deg)'
-      write(*,*) '(17)     peri_day : periastron date (sols since Ls=0)'
-      write(*,*) '(15)     periastr : min. star-planet dist (UA)'
-      write(*,*) '(16)     apoastr  : max. star-planet (UA)'
-      write(*,*) '(14)     year_day : length of year (in sols)'
-      write(*,*) '(5)      rad      : radius of the planet (m)'
-      write(*,*) '(6)      omeg     : planet rotation rate (rad/s)'
-      write(*,*) '(7)      g        : gravity (m/s2)'
-      write(*,*) '(8)      mugaz    : molecular mass '
-      write(*,*) '                       of the atmosphere (g/mol)'
-      write(*,*) '(9)      rcp      : r/Cp'
-      write(*,*) '(8)+(9)  calc_cpp_mugaz : r/Cp and mugaz '
-      write(*,*) '                 computed from gases.def'
-      write(*,*) '(10)     daysec   : length of a sol (s)'
-      write(*,*)
- 
- 
-      do while(modif(1:1).ne.'hello')
-        write(*,*)
-        write(*,*)
-        write(*,*) 'Changes to perform ?'
-        write(*,*) '   (enter keyword or return )'
-        write(*,*)
-        read(*,fmt='(a20)') modif
-        if (modif(1:1) .eq. ' ') goto 999
- 
-        write(*,*)
-        write(*,*) modif(1:len_trim(modif)) , ' : '
-
-        if (modif(1:len_trim(modif)) .eq. 'day_ini') then
-          write(*,*) 'current value:',day_ini
-          write(*,*) 'enter new value:'
- 101      read(*,*,iostat=ierr) day_ini
-          if(ierr.ne.0) goto 101
-          write(*,*) ' '
-          write(*,*) 'day_ini (new value):',day_ini
-
-        else if (modif(1:len_trim(modif)) .eq. 'z0') then
-          write(*,*) 'current value:',z0
-          write(*,*) 'enter new value:'
- 102      read(*,*,iostat=ierr) z0
-          if(ierr.ne.0) goto 102
-          write(*,*) ' '
-          write(*,*) ' z0 (new value):',z0
-
-        else if (modif(1:len_trim(modif)) .eq. 'emin_turb') then
-          write(*,*) 'current value:',emin_turb
-          write(*,*) 'enter new value:'
- 103      read(*,*,iostat=ierr) emin_turb
-          if(ierr.ne.0) goto 103
-          write(*,*) ' '
-          write(*,*) ' emin_turb (new value):',emin_turb
-
-        else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then
-          write(*,*) 'current value:',lmixmin
-          write(*,*) 'enter new value:'
- 104      read(*,*,iostat=ierr) lmixmin
-          if(ierr.ne.0) goto 104
-          write(*,*) ' '
-          write(*,*) ' lmixmin (new value):',lmixmin
-
-        else if (modif(1:len_trim(modif)) .eq. 'emissiv') then
-          write(*,*) 'current value:',emissiv
-          write(*,*) 'enter new value:'
- 105      read(*,*,iostat=ierr) emissiv
-          if(ierr.ne.0) goto 105
-          write(*,*) ' '
-          write(*,*) ' emissiv (new value):',emissiv
-
-        else if (modif(1:len_trim(modif)) .eq. 'emisice') then
-          write(*,*) 'current value emisice(1) North:',emisice(1)
-          write(*,*) 'enter new value:'
- 106      read(*,*,iostat=ierr) emisice(1)
-          if(ierr.ne.0) goto 106
-          write(*,*) 
-          write(*,*) ' emisice(1) (new value):',emisice(1)
-          write(*,*)
-
-          write(*,*) 'current value emisice(2) South:',emisice(2)
-          write(*,*) 'enter new value:'
- 107      read(*,*,iostat=ierr) emisice(2)
-          if(ierr.ne.0) goto 107
-          write(*,*) 
-          write(*,*) ' emisice(2) (new value):',emisice(2)
-
-        else if (modif(1:len_trim(modif)) .eq. 'iceradius') then
-          write(*,*) 'current value iceradius(1) North:',iceradius(1)
-          write(*,*) 'enter new value:'
- 110      read(*,*,iostat=ierr) iceradius(1)
-          if(ierr.ne.0) goto 110
-          write(*,*) 
-          write(*,*) ' iceradius(1) (new value):',iceradius(1)
-          write(*,*)
-
-          write(*,*) 'current value iceradius(2) South:',iceradius(2)
-          write(*,*) 'enter new value:'
- 111      read(*,*,iostat=ierr) iceradius(2)
-          if(ierr.ne.0) goto 111
-          write(*,*) 
-          write(*,*) ' iceradius(2) (new value):',iceradius(2)
-
-        else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then
-          write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
-          write(*,*) 'enter new value:'
- 112      read(*,*,iostat=ierr) dtemisice(1)
-          if(ierr.ne.0) goto 112
-          write(*,*) 
-          write(*,*) ' dtemisice(1) (new value):',dtemisice(1)
-          write(*,*)
-
-          write(*,*) 'current value dtemisice(2) South:',dtemisice(2)
-          write(*,*) 'enter new value:'
- 113      read(*,*,iostat=ierr) dtemisice(2)
-          if(ierr.ne.0) goto 113
-          write(*,*) 
-          write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
-
-        else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
-          write(*,*) 'current value:',obliquit
-          write(*,*) 'obliquit should be 25.19 on current Mars'
-          write(*,*) 'enter new value:'
- 115      read(*,*,iostat=ierr) obliquit
-          if(ierr.ne.0) goto 115
-          write(*,*) 
-          write(*,*) ' obliquit (new value):',obliquit
-
-        else if (modif(1:len_trim(modif)) .eq. 'peri_day') then
-          write(*,*) 'current value:',peri_day
-          write(*,*) 'peri_day should be 485 on current Mars'
-          write(*,*) 'enter new value:'
- 116      read(*,*,iostat=ierr) peri_day
-          if(ierr.ne.0) goto 116
-          write(*,*) 
-          write(*,*) ' peri_day (new value):',peri_day
-
-        else if (modif(1:len_trim(modif)) .eq. 'periastr') then
-          write(*,*) 'current value:',periastr
-          write(*,*) 'periastr should be 206.66 on present-day Mars'
-          write(*,*) 'enter new value:'
- 117      read(*,*,iostat=ierr) periastr
-          if(ierr.ne.0) goto 117
-          write(*,*) 
-          write(*,*) ' periastr (new value):',periastr
- 
-        else if (modif(1:len_trim(modif)) .eq. 'apoastr') then
-          write(*,*) 'current value:',apoastr
-          write(*,*) 'apoastr should be 249.22 on present-day Mars'
-          write(*,*) 'enter new value:'
- 118      read(*,*,iostat=ierr) apoastr
-          if(ierr.ne.0) goto 118
-          write(*,*) 
-          write(*,*) ' apoastr (new value):',apoastr
- 
-        else if (modif(1:len_trim(modif)) .eq. 'volcapa') then
-          write(*,*) 'current value:',volcapa
-          write(*,*) 'enter new value:'
- 119      read(*,*,iostat=ierr) volcapa
-          if(ierr.ne.0) goto 119
-          write(*,*) 
-          write(*,*) ' volcapa (new value):',volcapa
-        
-        else if (modif(1:len_trim(modif)).eq.'rad') then
-          write(*,*) 'current value:',rad
-          write(*,*) 'enter new value:'
- 120      read(*,*,iostat=ierr) rad
-          if(ierr.ne.0) goto 120
-          write(*,*) 
-          write(*,*) ' rad (new value):',rad
-
-        else if (modif(1:len_trim(modif)).eq.'omeg') then
-          write(*,*) 'current value:',omeg
-          write(*,*) 'enter new value:'
- 121      read(*,*,iostat=ierr) omeg
-          if(ierr.ne.0) goto 121
-          write(*,*) 
-          write(*,*) ' omeg (new value):',omeg
-        
-        else if (modif(1:len_trim(modif)).eq.'g') then
-          write(*,*) 'current value:',g
-          write(*,*) 'enter new value:'
- 122      read(*,*,iostat=ierr) g
-          if(ierr.ne.0) goto 122
-          write(*,*) 
-          write(*,*) ' g (new value):',g
-
-        else if (modif(1:len_trim(modif)).eq.'mugaz') then
-          write(*,*) 'current value:',mugaz
-          write(*,*) 'enter new value:'
- 123      read(*,*,iostat=ierr) mugaz
-          if(ierr.ne.0) goto 123
-          write(*,*) 
-          write(*,*) ' mugaz (new value):',mugaz
-          r=8.314511/(mugaz/1000.0)
-          write(*,*) ' R (new value):',r
-
-        else if (modif(1:len_trim(modif)).eq.'rcp') then
-          write(*,*) 'current value:',rcp
-          write(*,*) 'enter new value:'
- 124      read(*,*,iostat=ierr) rcp
-          if(ierr.ne.0) goto 124
-          write(*,*) 
-          write(*,*) ' rcp (new value):',rcp
-          r=8.314511/(mugaz/1000.0)
-          cpp=r/rcp
-          write(*,*) ' cpp (new value):',cpp
-
-        else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then
-          write(*,*) 'current value rcp, mugaz:',rcp,mugaz
-          check_cpp_match=.false.
-	  force_cpp=.false.
-	  call su_gases
-	  call calc_cpp_mugaz
-          write(*,*) 
-          write(*,*) ' cpp (new value):',cpp
-          write(*,*) ' mugaz (new value):',mugaz
-          r=8.314511/(mugaz/1000.0)
-          rcp=r/cpp
-          write(*,*) ' rcp (new value):',rcp
-	  
-        else if (modif(1:len_trim(modif)).eq.'daysec') then
-          write(*,*) 'current value:',daysec
-          write(*,*) 'enter new value:'
- 125      read(*,*,iostat=ierr) daysec
-          if(ierr.ne.0) goto 125
-          write(*,*) 
-          write(*,*) ' daysec (new value):',daysec
-
-!         added by RW!
-        else if (modif(1:len_trim(modif)).eq.'year_day') then
-          write(*,*) 'current value:',year_day
-          write(*,*) 'enter new value:' 
- 126      read(*,*,iostat=ierr) year_day
-          if(ierr.ne.0) goto 126
-          write(*,*)
-          write(*,*) ' year_day (new value):',year_day
-
-        endif
-      enddo ! of do while(modif(1:1).ne.'hello')
-
- 999  continue
-
-c-----------------------------------------------------------------------
-c	Write values of physical constants after modifications
-c-----------------------------------------------------------------------
- 
-      write(*,*) '*****************************************************'
-      write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes'
-      write(*,*) '*****************************************************'
-      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
-      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
-      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
-      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
-      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
-      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
-      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
-      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
-      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
-      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
- 
-      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
-      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
-      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
-      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
-      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
- 
-      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
-      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
-      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
- 
-      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
-      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
-      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
-      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
-      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
-      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
-      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
- 
-      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
-
-      write(*,*)  
-      write(*,*) 
-
-      ENDIF ! of if (Lmodif == 1)
-
-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
-
-
-      end
Index: trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90	(revision 1670)
+++ trunk/LMDZ.TITAN/libf/phytitan/tabfi_mod.F90	(revision 1670)
@@ -0,0 +1,585 @@
+MODULE tabfi_mod
+
+IMPLICIT NONE
+
+CONTAINS
+
+!=======================================================================
+      SUBROUTINE tabfi(ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad, &
+                       p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+!=======================================================================
+!
+!   C. Hourdin 15/11/96
+!
+!   Object:        Lecture du tab_cntrl physique dans un fichier 
+!   ------            et initialisation des constantes physiques
+!
+!   Arguments:
+!   ----------
+!
+!     Inputs:
+!     ------
+!
+!      - nid:    unitne logique du fichier ou on va lire le tab_cntrl    
+!                      (ouvert dans le programme appellant) 
+!
+!                 si nid=0:
+!                       pas de lecture du tab_cntrl mais
+!                       Valeurs par default des constantes physiques
+!        
+!      - tab0:    Offset de tab_cntrl a partir duquel sont ranges 
+!                  les parametres physiques (50 pour start_archive)
+!
+!      - Lmodif:  si on souhaite modifier les constantes  Lmodif = 1 = TRUE
+!
+!
+!     Outputs:
+!     --------
+!
+!      - day_ini: tab_cntrl(tab0+3) (Dans les cas ou l'on souhaite
+!                              comparer avec le day_ini dynamique)
+!
+!      - lmax:    tab_cntrl(tab0+2) (pour test avec nlayer)
+!
+!      - p_rad
+!      - p_omeg   !
+!      - p_g      ! Constantes physiques ayant des 
+!      - p_mugaz  ! homonymes dynamiques
+!      - p_daysec !
+!
+!=======================================================================
+! to use  'getin_p'
+      use ioipsl_getin_p_mod, only: getin_p
+
+      use surfdat_h, only: emisice, iceradius, dtemisice, &
+                           emissiv
+      use comsoil_h, only: volcapa
+      use iostart, only: get_var
+      use mod_phys_lmdz_para, only: is_parallel
+      use planete_mod, only: year_day, periastr, apoastr, peri_day, &
+                             obliquit, z0, lmixmin, emin_turb
+      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
+      use time_phylmdz_mod, only: dtphys, daysec
+      use callkeys_mod, only: check_cpp_match,force_cpp
+      implicit none
+ 
+      include "netcdf.inc"
+
+!-----------------------------------------------------------------------
+!   Declarations
+!-----------------------------------------------------------------------
+
+! Arguments
+! ---------
+      INTEGER,INTENT(IN) :: ngrid,nid,tab0
+      INTEGER*4,INTENT(OUT) :: day_ini
+      INTEGER,INTENT(IN) :: Lmodif
+      INTEGER,INTENT(OUT) :: lmax
+      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
+
+! Variables
+! ---------
+      INTEGER,PARAMETER :: length=100
+      REAL tab_cntrl(length) ! array in which are stored the run's parameters
+      INTEGER  ierr,nvarid
+      INTEGER size
+      CHARACTER modif*20
+      LOGICAL :: found
+      CHARACTER(len=5) :: modname="tabfi"
+      
+      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
+
+      IF (nid.eq.0) then
+!-----------------------------------------------------------------------
+!  Initialization of various physical constants to defaut values (nid = 0 case)
+!-----------------------------------------------------------------------
+        tab_cntrl(:)=0
+        lmax=0 ! not used anyways
+        !day_ini already set via inifis
+        time=0 
+! Informations about planet for dynamics and physics
+        ! rad,cpp,g,r,rcp already initialized by inifis
+        omeg=-999.
+        call getin_p("omega",omeg)
+        if (omeg.eq.-999.) then
+          call abort_physic(modname,"Missing value for omega in def files!",1)
+        endif
+        mugaz=(8.3144621/r)*1.E3
+        ! daysec already set by inifis
+        ! dtphys alread set by inifis
+! Informations about planet for the physics only
+        year_day=-999. ! length of year, in standard days
+        call getin_p("year_day",year_day)
+        if (year_day.eq.-999.) then
+          call abort_physic(modname, &
+               "Missing value for year_day in def files!",1)
+        endif
+        periastr=-999.
+        call getin_p("periastron",periastr)
+        if (periastr.eq.-999.) then
+          call abort_physic(modname, &
+               "Missing value for periastron in def files!",1)
+        endif
+        apoastr=-999.
+        call getin_p("apoastron",apoastr)
+        if (apoastr.eq.-999.) then
+          call abort_physic(modname, &
+               "Missing value for apoastron in def files!",1)
+        endif
+        peri_day=-999.
+        call getin_p("periastron_day",peri_day)
+        if (peri_day.eq.-999.) then
+          call abort_physic(modname, &
+               "Missing value for periastron date in def files!",1)
+        endif
+        obliquit=-999.
+        call getin_p("obliquity",obliquit)
+        if (obliquit.eq.-999.) then
+          call abort_physic(modname, &
+               "Missing value for obliquity in def files!",1)
+        endif
+! boundary layer and turbulence
+        z0=1.e-2 ! surface roughness length (m)
+        lmixmin=30
+        emin_turb=1.e-6
+! optical properties of polar caps and ground emissivity
+        emisice(:)=0
+        emissiv=0
+        iceradius(:)=1.e-6 ! mean scat radius of CO2 snow
+        dtemisice(:)=0 !time scale for snow metamorphism
+        volcapa=1000000 ! volumetric heat capacity of subsurface
+
+      ELSE
+!-----------------------------------------------------------------------
+!  Initialization of physical constants by reading array tab_cntrl(:)
+!		which contains these parameters	(nid != 0 case)
+!-----------------------------------------------------------------------
+! Read 'controle' array
+!
+
+       call get_var("controle",tab_cntrl,found)
+       if (.not.found) then
+         call abort_physic(modname,"Failed reading <controle> array",1)
+       else
+         write(*,*)'tabfi: tab_cntrl',tab_cntrl
+       endif
+!
+!  Initialization of some physical constants
+! informations on physics grid
+!      if(ngrid.ne.tab_cntrl(tab0+1)) then
+!         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
+!         print*,tab_cntrl(tab0+1),ngrid
+!      endif
+      lmax = nint(tab_cntrl(tab0+2))
+      day_ini = tab_cntrl(tab0+3)
+      time = tab_cntrl(tab0+4)
+      write (*,*) 'IN tabfi day_ini=',day_ini
+! Informations about planet for dynamics and physics
+      rad = tab_cntrl(tab0+5)
+      omeg = tab_cntrl(tab0+6)
+      g = tab_cntrl(tab0+7)
+      mugaz = tab_cntrl(tab0+8)
+      rcp = tab_cntrl(tab0+9)
+      cpp=(8.314511/(mugaz/1000.0))/rcp
+      daysec = tab_cntrl(tab0+10)
+      dtphys = tab_cntrl(tab0+11)
+! Informations about planet for the physics only
+      year_day = tab_cntrl(tab0+14)
+      periastr = tab_cntrl(tab0+15)
+      apoastr = tab_cntrl(tab0+16)
+      peri_day = tab_cntrl(tab0+17)
+      obliquit = tab_cntrl(tab0+18)
+! boundary layer and turbulence
+      z0 = tab_cntrl(tab0+19)
+      lmixmin = tab_cntrl(tab0+20)
+      emin_turb = tab_cntrl(tab0+21)
+! optical properties of polar caps and ground emissivity
+      emisice(1) = tab_cntrl(tab0+24)
+      emisice(2) = tab_cntrl(tab0+25)
+      emissiv    = tab_cntrl(tab0+26)
+      iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
+      iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
+      dtemisice(1)= tab_cntrl(tab0+33) !time scale for snow metamorphism (north)
+      dtemisice(2)= tab_cntrl(tab0+34) !time scale for snow metamorphism (south)
+! soil properties
+      volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
+!-----------------------------------------------------------------------
+!	Save some constants for later use (as routine arguments)
+!-----------------------------------------------------------------------
+      p_omeg = omeg
+      p_g = g
+      p_cpp = cpp
+      p_mugaz = mugaz
+      p_daysec = daysec
+      p_rad=rad
+
+      ENDIF    ! end of (nid = 0) 
+
+!-----------------------------------------------------------------------
+!	Write physical constants to output before modifying them
+!-----------------------------------------------------------------------
+ 
+   6  FORMAT(a20,e15.6,e15.6)
+   5  FORMAT(a20,f12.2,f12.2)
+ 
+      write(*,*) '*****************************************************'
+      write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
+      write(*,*) '*****************************************************'
+      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
+      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
+      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
+      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
+      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
+      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
+      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
+      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
+      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
+      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
+
+      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
+      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
+      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
+      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
+      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
+
+      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
+      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
+      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
+
+      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
+      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
+      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
+      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
+      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
+      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
+      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
+
+      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
+
+      write(*,*)
+      write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif
+
+!-----------------------------------------------------------------------
+!	 Modifications...
+! NB: Modifying controls should only be done by newstart, and in seq mode
+      if ((Lmodif.eq.1).and.is_parallel) then
+        write(*,*) "tabfi: Error modifying tab_control should", &
+                   " only happen in serial mode (eg: by newstart)"
+        stop
+      endif
+!-----------------------------------------------------------------------
+
+      IF(Lmodif.eq.1) then
+
+      write(*,*)
+      write(*,*) 'Change values in tab_cntrl ? :'
+      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+      write(*,*) '(Current values given above)'
+      write(*,*)
+      write(*,*) '(3)          day_ini : Initial day (=0 at Ls=0)'
+      write(*,*) '(19)              z0 :  surface roughness (m)'
+      write(*,*) '(21)       emin_turb :  minimal energy (PBL)'
+      write(*,*) '(20)         lmixmin : mixing length (PBL)'
+      write(*,*) '(26)         emissiv : ground emissivity'
+      write(*,*) '(24 et 25)   emisice : CO2 ice max emissivity '
+      write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
+      write(*,*) '(33 et 34) dtemisice : time scale for snow metamorphism'
+      write(*,*) '(35)      volcapa : soil volumetric heat capacity'
+      write(*,*) '(18)     obliquit : planet obliquity (deg)'
+      write(*,*) '(17)     peri_day : periastron date (sols since Ls=0)'
+      write(*,*) '(15)     periastr : min. star-planet dist (UA)'
+      write(*,*) '(16)     apoastr  : max. star-planet (UA)'
+      write(*,*) '(14)     year_day : length of year (in sols)'
+      write(*,*) '(5)      rad      : radius of the planet (m)'
+      write(*,*) '(6)      omeg     : planet rotation rate (rad/s)'
+      write(*,*) '(7)      g        : gravity (m/s2)'
+      write(*,*) '(8)      mugaz    : molecular mass '
+      write(*,*) '                       of the atmosphere (g/mol)'
+      write(*,*) '(9)      rcp      : r/Cp'
+      write(*,*) '(8)+(9)  calc_cpp_mugaz : r/Cp and mugaz '
+      write(*,*) '                 computed from gases.def'
+      write(*,*) '(10)     daysec   : length of a sol (s)'
+      write(*,*)
+ 
+ 
+      do while(modif(1:1).ne.'hello')
+        write(*,*)
+        write(*,*)
+        write(*,*) 'Changes to perform ?'
+        write(*,*) '   (enter keyword or return )'
+        write(*,*)
+        read(*,fmt='(a20)') modif
+        if (modif(1:1) .eq. ' ') goto 999
+ 
+        write(*,*)
+        write(*,*) modif(1:len_trim(modif)) , ' : '
+
+        if (modif(1:len_trim(modif)) .eq. 'day_ini') then
+          write(*,*) 'current value:',day_ini
+          write(*,*) 'enter new value:'
+ 101      read(*,*,iostat=ierr) day_ini
+          if(ierr.ne.0) goto 101
+          write(*,*) ' '
+          write(*,*) 'day_ini (new value):',day_ini
+
+        else if (modif(1:len_trim(modif)) .eq. 'z0') then
+          write(*,*) 'current value:',z0
+          write(*,*) 'enter new value:'
+ 102      read(*,*,iostat=ierr) z0
+          if(ierr.ne.0) goto 102
+          write(*,*) ' '
+          write(*,*) ' z0 (new value):',z0
+
+        else if (modif(1:len_trim(modif)) .eq. 'emin_turb') then
+          write(*,*) 'current value:',emin_turb
+          write(*,*) 'enter new value:'
+ 103      read(*,*,iostat=ierr) emin_turb
+          if(ierr.ne.0) goto 103
+          write(*,*) ' '
+          write(*,*) ' emin_turb (new value):',emin_turb
+
+        else if (modif(1:len_trim(modif)) .eq. 'lmixmin') then
+          write(*,*) 'current value:',lmixmin
+          write(*,*) 'enter new value:'
+ 104      read(*,*,iostat=ierr) lmixmin
+          if(ierr.ne.0) goto 104
+          write(*,*) ' '
+          write(*,*) ' lmixmin (new value):',lmixmin
+
+        else if (modif(1:len_trim(modif)) .eq. 'emissiv') then
+          write(*,*) 'current value:',emissiv
+          write(*,*) 'enter new value:'
+ 105      read(*,*,iostat=ierr) emissiv
+          if(ierr.ne.0) goto 105
+          write(*,*) ' '
+          write(*,*) ' emissiv (new value):',emissiv
+
+        else if (modif(1:len_trim(modif)) .eq. 'emisice') then
+          write(*,*) 'current value emisice(1) North:',emisice(1)
+          write(*,*) 'enter new value:'
+ 106      read(*,*,iostat=ierr) emisice(1)
+          if(ierr.ne.0) goto 106
+          write(*,*) 
+          write(*,*) ' emisice(1) (new value):',emisice(1)
+          write(*,*)
+
+          write(*,*) 'current value emisice(2) South:',emisice(2)
+          write(*,*) 'enter new value:'
+ 107      read(*,*,iostat=ierr) emisice(2)
+          if(ierr.ne.0) goto 107
+          write(*,*) 
+          write(*,*) ' emisice(2) (new value):',emisice(2)
+
+        else if (modif(1:len_trim(modif)) .eq. 'iceradius') then
+          write(*,*) 'current value iceradius(1) North:',iceradius(1)
+          write(*,*) 'enter new value:'
+ 110      read(*,*,iostat=ierr) iceradius(1)
+          if(ierr.ne.0) goto 110
+          write(*,*) 
+          write(*,*) ' iceradius(1) (new value):',iceradius(1)
+          write(*,*)
+
+          write(*,*) 'current value iceradius(2) South:',iceradius(2)
+          write(*,*) 'enter new value:'
+ 111      read(*,*,iostat=ierr) iceradius(2)
+          if(ierr.ne.0) goto 111
+          write(*,*) 
+          write(*,*) ' iceradius(2) (new value):',iceradius(2)
+
+        else if (modif(1:len_trim(modif)) .eq. 'dtemisice') then
+          write(*,*) 'current value dtemisice(1) North:',dtemisice(1)
+          write(*,*) 'enter new value:'
+ 112      read(*,*,iostat=ierr) dtemisice(1)
+          if(ierr.ne.0) goto 112
+          write(*,*) 
+          write(*,*) ' dtemisice(1) (new value):',dtemisice(1)
+          write(*,*)
+
+          write(*,*) 'current value dtemisice(2) South:',dtemisice(2)
+          write(*,*) 'enter new value:'
+ 113      read(*,*,iostat=ierr) dtemisice(2)
+          if(ierr.ne.0) goto 113
+          write(*,*) 
+          write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
+
+        else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
+          write(*,*) 'current value:',obliquit
+          write(*,*) 'obliquit should be 25.19 on current Mars'
+          write(*,*) 'enter new value:'
+ 115      read(*,*,iostat=ierr) obliquit
+          if(ierr.ne.0) goto 115
+          write(*,*) 
+          write(*,*) ' obliquit (new value):',obliquit
+
+        else if (modif(1:len_trim(modif)) .eq. 'peri_day') then
+          write(*,*) 'current value:',peri_day
+          write(*,*) 'peri_day should be 485 on current Mars'
+          write(*,*) 'enter new value:'
+ 116      read(*,*,iostat=ierr) peri_day
+          if(ierr.ne.0) goto 116
+          write(*,*) 
+          write(*,*) ' peri_day (new value):',peri_day
+
+        else if (modif(1:len_trim(modif)) .eq. 'periastr') then
+          write(*,*) 'current value:',periastr
+          write(*,*) 'periastr should be 206.66 on present-day Mars'
+          write(*,*) 'enter new value:'
+ 117      read(*,*,iostat=ierr) periastr
+          if(ierr.ne.0) goto 117
+          write(*,*) 
+          write(*,*) ' periastr (new value):',periastr
+ 
+        else if (modif(1:len_trim(modif)) .eq. 'apoastr') then
+          write(*,*) 'current value:',apoastr
+          write(*,*) 'apoastr should be 249.22 on present-day Mars'
+          write(*,*) 'enter new value:'
+ 118      read(*,*,iostat=ierr) apoastr
+          if(ierr.ne.0) goto 118
+          write(*,*) 
+          write(*,*) ' apoastr (new value):',apoastr
+ 
+        else if (modif(1:len_trim(modif)) .eq. 'volcapa') then
+          write(*,*) 'current value:',volcapa
+          write(*,*) 'enter new value:'
+ 119      read(*,*,iostat=ierr) volcapa
+          if(ierr.ne.0) goto 119
+          write(*,*) 
+          write(*,*) ' volcapa (new value):',volcapa
+        
+        else if (modif(1:len_trim(modif)).eq.'rad') then
+          write(*,*) 'current value:',rad
+          write(*,*) 'enter new value:'
+ 120      read(*,*,iostat=ierr) rad
+          if(ierr.ne.0) goto 120
+          write(*,*) 
+          write(*,*) ' rad (new value):',rad
+
+        else if (modif(1:len_trim(modif)).eq.'omeg') then
+          write(*,*) 'current value:',omeg
+          write(*,*) 'enter new value:'
+ 121      read(*,*,iostat=ierr) omeg
+          if(ierr.ne.0) goto 121
+          write(*,*) 
+          write(*,*) ' omeg (new value):',omeg
+        
+        else if (modif(1:len_trim(modif)).eq.'g') then
+          write(*,*) 'current value:',g
+          write(*,*) 'enter new value:'
+ 122      read(*,*,iostat=ierr) g
+          if(ierr.ne.0) goto 122
+          write(*,*) 
+          write(*,*) ' g (new value):',g
+
+        else if (modif(1:len_trim(modif)).eq.'mugaz') then
+          write(*,*) 'current value:',mugaz
+          write(*,*) 'enter new value:'
+ 123      read(*,*,iostat=ierr) mugaz
+          if(ierr.ne.0) goto 123
+          write(*,*) 
+          write(*,*) ' mugaz (new value):',mugaz
+          r=8.314511/(mugaz/1000.0)
+          write(*,*) ' R (new value):',r
+
+        else if (modif(1:len_trim(modif)).eq.'rcp') then
+          write(*,*) 'current value:',rcp
+          write(*,*) 'enter new value:'
+ 124      read(*,*,iostat=ierr) rcp
+          if(ierr.ne.0) goto 124
+          write(*,*) 
+          write(*,*) ' rcp (new value):',rcp
+          r=8.314511/(mugaz/1000.0)
+          cpp=r/rcp
+          write(*,*) ' cpp (new value):',cpp
+
+        else if (modif(1:len_trim(modif)).eq.'calc_cpp_mugaz') then
+          write(*,*) 'current value rcp, mugaz:',rcp,mugaz
+          check_cpp_match=.false.
+	  force_cpp=.false.
+	  call su_gases
+	  call calc_cpp_mugaz
+          write(*,*) 
+          write(*,*) ' cpp (new value):',cpp
+          write(*,*) ' mugaz (new value):',mugaz
+          r=8.314511/(mugaz/1000.0)
+          rcp=r/cpp
+          write(*,*) ' rcp (new value):',rcp
+	  
+        else if (modif(1:len_trim(modif)).eq.'daysec') then
+          write(*,*) 'current value:',daysec
+          write(*,*) 'enter new value:'
+ 125      read(*,*,iostat=ierr) daysec
+          if(ierr.ne.0) goto 125
+          write(*,*) 
+          write(*,*) ' daysec (new value):',daysec
+
+!         added by RW!
+        else if (modif(1:len_trim(modif)).eq.'year_day') then
+          write(*,*) 'current value:',year_day
+          write(*,*) 'enter new value:' 
+ 126      read(*,*,iostat=ierr) year_day
+          if(ierr.ne.0) goto 126
+          write(*,*)
+          write(*,*) ' year_day (new value):',year_day
+
+        endif
+      enddo ! of do while(modif(1:1).ne.'hello')
+
+ 999  continue
+
+!----------------------------------------------------------------------
+!	Write values of physical constants after modifications
+!-----------------------------------------------------------------------
+ 
+      write(*,*) '*****************************************************'
+      write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes'
+      write(*,*) '*****************************************************'
+      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
+      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
+      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
+      write(*,5) '(5)             rad',tab_cntrl(tab0+5),rad
+      write(*,5) '(10)         daysec',tab_cntrl(tab0+10),daysec
+      write(*,6) '(6)            omeg',tab_cntrl(tab0+6),omeg
+      write(*,5) '(7)               g',tab_cntrl(tab0+7),g
+      write(*,5) '(8)           mugaz',tab_cntrl(tab0+8),mugaz
+      write(*,5) '(9)             rcp',tab_cntrl(tab0+9),rcp
+      write(*,6) '(11)        dtphys?',tab_cntrl(tab0+11),dtphys
+ 
+      write(*,5) '(14)       year_day',tab_cntrl(tab0+14),year_day
+      write(*,5) '(15)       periastr',tab_cntrl(tab0+15),periastr
+      write(*,5) '(16)        apoastr',tab_cntrl(tab0+16),apoastr
+      write(*,5) '(17)       peri_day',tab_cntrl(tab0+17),peri_day
+      write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
+ 
+      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
+      write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
+      write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
+ 
+      write(*,5) '(26)        emissiv',tab_cntrl(tab0+26),emissiv
+      write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
+      write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
+      write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
+      write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
+      write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
+      write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
+ 
+      write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
+
+      write(*,*)  
+      write(*,*) 
+
+      ENDIF ! of if (Lmodif == 1)
+
+!-----------------------------------------------------------------------
+!	Save some constants for later use (as routine arguments)
+!-----------------------------------------------------------------------
+      p_omeg = omeg
+      p_g = g
+      p_cpp = cpp
+      p_mugaz = mugaz
+      p_daysec = daysec
+      p_rad=rad
+
+
+      END SUBROUTINE tabfi
+
+end module tabfi_mod
